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Abstract 


The  kinematics  and  flow  field  development  of  two  wing  sections  undergoing  stall  flutter  limit 
cycle  oscillations  was  thoroughly  investigated  with  the  goal  of  elucidating  future  pathways 
for  successfully  implementing  aerodynamic  flow  control  on  flexible  wing  structures.  A  cyber¬ 
physical  system  was  implemented  with  each  wing  section  to  prescribe  and  tune,  in-situ,  the 
torsional  dynamics  of  the  wing  structure.  The  first  wing  section  studied  was  a  flexible,  semi- 
span,  rectangular  wing  section  composed  of  a  NACA  0018  profile  and  having  a  semi-span 
aspect  ratio  of  six.  When  tuned  appropriately,  the  wing  was  able  to  maintain  repeatable 
stall  flutter  limit  cycle  oscillations,  which  were  dominated  by  primarily  torsional  motions. 
The  torsional  and  bending  kinematics  were  found  to  be  accurately  represented  the  first 
three  eigenmodes  of  a  fixed-free  beam  in  free  vibration.  When  analyzing  the  flowfield  it  was 
found  that  a  coherent  dynamic  stall  vortex  was  not  produced  in  contrast  to  prior  work  on 
rigid  wings.  Instead  the  wing  stall  and  flow  separation  was  dominated  by  the  growth  and 
collapse  of  a  small  parabolic  region  of  flow  separation  centered  around  75%  span  location 
and  at  the  wing  trailing  edge.  The  second  wing  section  studied  was  a  rigid,  full-span, 
rectangular  wing  section  composed  of  a  NACA  0018  profile  and  having  an  aspect  ration 
of  4.  Instead  of  undergoing  a  twisting  deformation  this  wing  was  constrained  to  a  rigid- 
body  pitching  motion  in  the  wind  tunnel  and  was  used  to  identify  the  differences  between 
prescribed  sinusoidal  pitching  motion  and  passively  sustained  stall  flutter.  Specifically, 
it  was  found  that  stall  flutter  displays  a  slower  pitch-up  rate  as  compared  to  classically 
prescribed  sinusoidal  motions.  This  slower  pitch-up  rate  reduced  the  maximum  magnitude 
of  the  induced  pitching  moment  and  the  coherence  of  the  primary  dynamic  stall  vortex  that 
was  formed. 
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1  Summary 

The  following  final  performance  report  documents  the  technical  research  conducted  by  the 
principle  investigator  (PI)  and  his  graduate  research  assistant  (GRA),  Ethan  Culler,  during 
the  full  performance  period  (September  24,  2014  through  May  31,  2017)  of  cooperative 
agreement  number  FA7000-14-2-0018  entitled  Flow  Control  of  Flexible  Structures.  This 
report  focuses  on  the  technical  progress  and  achievements  under  this  agreement,  with  the 
financial  progress  reported  separately  in  quarterly  financial  reports  as  documented  in  the 
cooperative  agreement. 

The  final  report  is  structured  to  address  the  technical  objectives  that  were  laid  out  in 
the  original  proposal.  As  a  result  the  original  scientific  objectives  are  summarized  here  in 
Section  2.1  which  is  followed  by  a  brief  discussion  of  deviations  from  the  proposed  effort  in 
Section  2.2.  It  should  be  noted  that  an  extensive  literature  review  is  not  included  in  this 
report  for  brevity;  if  readers  are  interested  in  a  more  detailed  discussion  of  the  prior  work 
they  are  referred  to  journal  and  conference  publications  that  were  published  as  a  result 
of  this  effort  which  are  listed  in  Section  7.  Instead  this  report  will  focus  on  the  scientific 
results  and  accomplishments  made  under  this  research  effort.  As  a  result  the  experimental 
methodology  is  discussed  in  detail  in  Section  3  which  is  followed  by  a  thorough  discussion 
of  the  experimental  results  and  findings  in  Section  4.  Additionally,  the  results  section  is 
divided  into  to  subsections  which  align  with  the  two  major  phases  of  the  proposed  effort. 
The  first  phase,  (Section  4.1)  focused  on  the  detailed  analysis  of  the  wing  kinematics  and 
flowfield  development  for  the  flexible  cyber-physical  wing  system  undergoing  stall  flutter. 
These  experiments  were  carried  out  in  the  Subsonic  Wind  Tunnel  at  the  United  States  Air 
Force  Academy  (USAFA)  and  the  results  were  then  analyzed  by  the  PI  and  GRA  back  at  the 
University  of  Colorado  Boulder  (UCB).  The  second  phase  of  the  research  effort  (Section  4.2) 
focused  on  elucidating  the  differences  between  a  rigid  wing  undergoing  sinusoidal  prescribed 
motion  and  passivily  sustained  stall  flutter  limit  cycle  oscillations  (LCOs).  These  experi¬ 
ments  were  carried  out  at  both  the  USAFA  and  UCB  and  again  the  results  were  analyzed 
by  the  PI  and  GRA  at  UCB.  The  detailed  discussion  of  the  results,  is  followed  by  a  brief 
summary  of  the  accomplishments  of  this  research  effort  (i.e.  publications,  presentations, 
and  awards)  in  Section  7.  The  a  summary  of  the  scientific  conclusions  from  the  research 
effort  is  made  in  Section  5,  which  is  followed  by  a  brief  discussion  of  recomendations  for 
future  work  in  Section  6. 

2  Introduction 

Fluid-structure  interactions  have  long  been  a  research  topic  of  great  interest.  However,  to 
date,  most  investigations  regard  the  fluid  dynamics  problem  as  too  complex  to  be  fully 
understood  and  therefore  resort  to  various  levels  of  modeling  [Dowell  and  Hall,  2001].  Fur¬ 
thermore,  even  for  the  most  complex  models,  the  main  research  goal  has  been  a  mathe¬ 
matical  description  of  the  structural  dynamics  with  the  flow  merely  providing  the  “forcing 
term”  [Barone  and  Payne,  2005].  With  recent  advances  in  computational  and  experimental 
technology,  however,  it  is  now  possible  to  consider  what  can  be  learned  about  the  unsteady 
physical  processes  in  the  flowfield  that  eventually  result  in  structural  deformations.  Ulti¬ 
mately,  research  in  fluid-structure  interaction  can  show  ways  to  tailor  the  flowfield  (using 
feedback  flow  control)  or  the  structure  (using  adaptive  materials  and  active  structural  con- 
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trol)  to  accomplish  a  prescribed  objective. 

A  comparison  of  the  dynamic  behavior  of  fluid  flowing  over  a  flexible  structure  and  the 
structure  itself  reveals  that  there  are  fundamentally  different  time  scales  in  the  fluid  and 
structural  dynamics,  where  the  fluid  dynamic  processes  occur  on  a  significantly  faster  time 
scale  and  a  concomitantly  shorter  length  scale.  In  fact,  recent  research  unequivocally  shows 
that  the  response  of  a  structure  to  aerodynamic  loading  has  to  be  understood  as  the  integral 
effect  of  the  faster  and  smaller  scale  fluid  dynamics  occurring  in  the  flow  around  the  struc¬ 
ture  [Mangalam  et  al.,  2008].  Therefore,  sensing  and  controlling  the  fluid  dynamics  instead 
of  solely  measuring  and  controlling  the  structural  deformation  (e.g.  with  accelerometers) 
provides  two  key  advantages:  i)  early  detection  of  flowfield  changes  that  will  eventually 
result  in  structural  loads,  and  ii)  the  possibility  for  localized,  discrete  actuation  to  coun¬ 
teract  detrimental  flow  developments  before  they  result  in  significant  structural  loads  and 
deformations.  If  such  sensing  and  control  is  possible,  significant  improvements  to  engineer¬ 
ing  systems  can  be  expected.  Applications  such  as  renewable  energy  systems  (e.g.  wind 
turbines  or  ocean  energy  devices),  air  vehicle  aerodynamics  and  engines,  or  even  medical 
flows  (blood  flow,  respiration)  could  be  better  understood  and  therefore  their  performance 
could  be  improved. 

Aerodynamically,  fluid-structure  interactions  can  be  primarily  placed  in  two  categories: 
those  that  are  induced  by  lift  and  those  induced  by  drag.  For  lift-induced  fluid-structure 
interaction,  controlling  the  lift  distribution  provides  a  means  to  control  the  structural  defor¬ 
mation.  This  is  the  classical  aeroelasticity  problem,  for  which  passive  and  active  structural 
control  strategies  have  been  thoroughly  investigated  [Fung,  1955,  Librescu  and  Marzocca, 
2005].  As  a  result  this  represents  the  primary  focus  of  the  current  research. 

From  an  active  flow  control  perspective,  nature  provides  insight  into  how  aeroelastic  sys¬ 
tems  can  be  controlled:  with  sophisticated  methods  integrating  instantaneous  distributed 
flow  state  sensing  and  control  actuation.  Birds,  for  example,  use  their  feathers  as  both 
sensors  to  monitor  the  local  flow  state  and  as  actuators  to  achieve  local  flow  control.  Bats, 
on  the  other  hand,  use  fine  hairs  sparsely  distributed  on  their  wings  to  sense  the  air  flow 
and  morph  their  wings  through  altering  the  wing  geometry  and  the  tension  in  the  wing 
membrane,  as  necessary  to  achieve  a  desired  flight  path  in  the  instantaneous  prevailing  flow 
conditions  [Sterbing-D ’Angelo  et  al.,  2011].  Clearly,  nature  has  developed  tightly  coupled 
control  systems,  employing  flow  control  and  structural  control  simultaneously.  To  mimic 
these  tightly  coupled  feedback  control  systems,  sensors,  actuators,  and  control  strategies 
need  to  be  developed  to  achieve  control  of  flexible  structures  while  maintaining  aerody¬ 
namic  system  performance.  A  thorough  understanding  of  the  underlying  flowfield,  including 
the  effects  of  structural  flexibility,  is  crucially  important  for  the  eventual  development  of 
feedback  flow  control  on  flexible  structures  in  technical  applications. 

2.1  Proposed  Research  Effort 

A  collaborative  research  program  was  originally  proposed  with  the  goal  of  achieving  active 
flow  control  on  flexible  structures  through  a  progressive  sequence  of  research  objectives. 
That  said  the  proposed  effort  focused  on  understanding  the  fundamental  fluid  dynamics 
associated  with  the  fluid-structure  interaction  problem  with  the  goal  of  determining  how 
the  implementation  of  open-loop  fluidic  actuation  directly  alters  the  fluid  dynamics  and 
indirectly  alters  the  structural  dynamics.  To  achieve  this  goal  the  original  proposal  laid  out 
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the  following  scientific  objectives: 

1.  To  determine  if  the  same  three-dimensional  unsteady  vortical  structures  are  formed 
and  shed  from  a  passively  driven  (driven  solely  by  the  oncoming  flow)  aeroelastic  wing 
as  compared  to  an  actively  oscillated  wing  model  at  a  moderate  Reynolds  number. 

2.  To  investigate  the  role  of  the  flow  velocity  and  material  stiffness  in  the  formation  and 
growth  of  the  unsteady  structures  in  the  flowfield  and  thus  the  structural  oscillations. 

3.  To  determine  from  the  passive  (unforced)  characteristics  the  optimal  locations  for 
actuation  and  sensing  devices  on  the  flexible  wing  surface. 

4.  To  investigate  whether  fluidic  actuators  can  be  utilized  to  both  amplify  and  cancel 
the  formation  and  growth  of  the  unsteady  vortical  structures. 

5.  To  determine  if  the  different  vortex  wake  modes  (and  thus  aerodynamic  loads)  can  be 
generated  through  exciting  the  flowfield  and  structure  with  small  local  alterations  to 
the  unsteady  flowfield,  without  changing  the  free-streanr  velocity,  material  properties, 
or  mechanically  driving  the  wing. 

2.2  Deviations  from  Proposed  Effort 

As  highlighted  above,  the  proposed  research  represented  a  collaborative  effort  which  was 
closely  coordinated  with  research  efforts  in  the  Aeronautics  Laboratory  at  USAFA  lead  by 
Dr.  Jurgen  Seidel  and  Dr.  Casey  Fagley.  The  progress  of  the  collaborative  research  effort 
was  regularly  reviewed  throughout  the  performance  period  by  the  group  (i.e.  Dr.  Seidel, 
Dr.  Fagley,  the  GRA  and  PI).  The  short  term  plans  for  the  future  research  were  then  made 
by  the  group  based  upon  the  achieved  progress  and  resources  available  at  each  of  the  sites 
(i.e.  USAFA  and  UCB).  This  coordination  and  planning  allowed  the  research  effort  to  stay 
on  task,  but  also  allowed  for  deviations  for  the  proposed  effort  to  be  coordinated  as  needed. 

The  most  significant  deviation  from  the  proposed  effort  was  related  to  the  implementa¬ 
tion  of  flow  control  on  the  flexible  wing  system.  In  the  original  proposal  a  significant  amount 
of  effort  was  focused  on  implementing  flow  control  in  the  flexible  wing.  This  included  objec¬ 
tives  3,  4,  and  5  in  Section  2.1.  It  was  realized  early  on  by  the  group  that  accomplishing  all 
of  these  tasks  was  overly  ambitious  with  the  limited  availability  of  the  Subsonic  Wind  Tun¬ 
nel  (at  USAFA)  during  the  summer  testing  periods.  As  a  result  the,  GRA  and  PI  focused 
primarily  on  the  detailed  investigations  into  the  baseline  kinematics  and  flowfield  surround¬ 
ing  the  flexible  wing  in  stall  flutter  during  the  first  summer  testing  period  (May  -  August 
2015).  During  the  second  summer  testing  period  (May  -  August  2016),  the  Subsonic  Wind 
Tunnel  at  USAFA  was  being  upgraded  which  significantly  limited  the  available  testing  time 
for  the  proposed  research.  As  a  result  the  GRA  and  PI  worked  during  the  second  summer 
period  to  construct  a  second  cyber-physical  testing  assembly  at  UCB.  Once  operational, 
this  allowed  tests  to  be  carried  out  at  UCB  in  the  low-speed  research  wind  tunnel,  which 
had  less  limitations  on  availability.  Specifically,  the  testing  at  UCB  focused  on  quantifying 
the  differences  in  the  kinematics  and  flowfield  development  between  a  rigid  wing  body  in 
passively  sustained  stall  flutter  limit  cycle  oscillations  and  the  same  wing  actively  pitched 
(or  driven)  in  similar  periodic  oscillations.  Results  and  analysis  of  both  the  flexible  wing 
and  the  rigid  wing  problems  will  be  presented  in  detail  below  below  in  Section  4  after  a 
thorough  discussion  of  the  experimental  facilities  and  measurement  techniques  in  Section  3. 
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3  Methodology 

The  research  carried  out  as  a  part  of  this  effort  focused  entirely  on  experimental  measure¬ 
ments  of  two  distinctly  different  fluttering  wing  models,  namely:  1)  a  semi-span,  flexible 
wing  model  undergoing  spanwise  varying  twisting  oscillations  and  2)  a  full-span,  rigid  wing 
model  undergoing  spanwise  uniform  pitching  oscillations.  Note  that  a  cyber-physical  ap¬ 
proach  was  utilized  to  dynamically  prescribe  the  motions  for  both  of  these  model  config¬ 
urations.  Finally,  as  part  of  the  research  program  measurements  were  collected  in  both 
the  subsonic  wind  tunnel  at  USAFA  and  the  low-speed  research  wind  tunnel  at  UCB.  The 
following  subsections  provide  additional  details  about  both  of  these  facilities,  each  of  the 
model  configurations,  and  the  measurement  techniques  utilized. 

3.1  Subsonic  Wind  Tunnel,  USAFA 

The  United  States  Air  Force  Academy’s  Subsonic  Wind  Tunnel  is  located  in  the  Aeronautics 
Research  Center  and  was  fabricated  by  FluiDyne  Engineering.  It  is  a  single  leg  recirculating 
type  wind  tunnel  and  has  a  0.91  m  by  0.91  m  by  1.83  m  test  section.  It  is  capable  of  attaining 
speeds  up  to  Mach  numbers  of  M  =  0.5,  however  for  most  of  the  tests  discussed  here  the 
wind  speed  was  limited  to  U, ^  =  26  m/s  or  M  =  0.08.  The  tunnel  facility  is  instrumented 
with  a  MKS  Baratron  220C  transducer  to  measure  the  dynamic  pressure  in  the  test  section. 

3.2  Low-Speed  Research  Wind  Tunnel,  UCB 

The  Low-Speed  Research  Wind  Tunnel  in  the  Experimental  Aerodynamics  Laboratory  at 
the  University  of  Colorado  Boulder  was  designed  and  fabricated  by  Aerolab  LLC.  It  is  a 
open-return  facility  with  an  upstream  centrifugal  blower  and  has  a  test  section  of  0.76  m 
by  0.76  m  by  3.5  nr  in  length.  It  is  capable  of  attaining  speeds  up  to  Mach  numbers  of 
M  =  0.2,  however  again  for  most  of  the  tests  discussed  here  the  wind  speed  was  limited 
to  Coo  =  15  m/s  or  M  =  0.04.  The  tunnel  facility  is  instrumented  with  a  MKS  Baratron 
220D  transducer  to  measure  the  dynamic  pressure,  a  MKS  Baratron  220A  to  measure  the 
static  pressure,  and  a  T-Type  thermocouple  to  measure  the  total  temperature  within  the 
test  section. 

3.3  Cyber-Physical  System 

The  stall  flutter  limit  cycle  oscillations  of  the  two  wing  models  were  tuned  using  the  position- 
feedback  cyber-physical  control  system  originally  developed  by  Fagley  et  al.  [2015,  2016] 
and  pictured  in  Figure  1  for  the  flexible  wing  model.  The  rigid  wing  cyber-physical  system 
assembly  was  very  similar  and  is  pictured  in  a  horizontal  orientation  in  Figure  6.  For  the 
flexible  wing  experiments  conducted  at  USAFA,  the  actual  cyber-physical  system  developed 
and  employed  by  Fagley  et  al.  [2015,  2016]  was  utilized.  As  part  of  the  current  research 
effort,  a  replica  of  this  cyber-physical  system  was  built  at  UCB  and  utilized  for  the  rigid 
wing  experiments.  In  both  cases,  an  embedded  controller  was  designed  to  take  the  form  of 
a  second-order  mass-spring-damper  system  as  shown  in  Equation  (1).  The  angle  of  attack 
was  determined  from  a  16,  000  count  optical  encoder,  while  the  angular  velocity  and  angular 
acceleration  were  calculated  with  the  embedded  FPGA  on  a  National  Instruments  CRIO- 
9063.  This  produced  an  angle  of  attack  measurement  precision  of  0.02°  with  an  FPGA 
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Figure  1:  Illustration  of  the  cyber-physical  system  which  uses  feedback  control  to  simulate 
the  torsional  structural  dynamics  of  the  flexible  wing  model.  [Fagley  et  al.,  2015]. 


output  update  every  2  ms. 


JyO  +  riv0  +  kgd  =  TMotor-  (1) 

Using  this  cyber-physical  system  approach,  the  rotational  stiffness,  kg,  virtual  rotational 
damping  rjy ,  and  virtual  rotational  inertia,  Jy,  were  tuned  through  a  respective  change  in 
the  control  system  gains.  In  Equation  (1),  6  is  the  tip-twist  angle  (flexible  wing)  or  rigid- 
body  rotation  angle  (rigid  wing)  with  respect  to  its  unforced  or  wind  off  position  and  TMotor 
is  the  applied  torque.  Furthermore,  the  reference  position  for  9  is  the  geometric  angle  of 
attack,  aQ.  Thus  the  model’s  structural  dynamics  could  be  easily  adjusted  in-situ  without 
requiring  physical  modifications  to  the  wing  model.  The  second  order  dynamical  system  is 
prescribed  as  shown  by  the  control  system  block  diagram  in  Figure  2. 

3.4  Flexible  Wing  Model 

The  test  geometry  for  the  flexible  wing  model  consisted  of  a  finite  span  NACA  0018  wing 
model  cantilevered  up  from  a  splitter  plate  mounted  to  the  wind  tunnel  floor  as  shown  in 
Figure  1.  Note  that  the  splitter  plate  was  utilized  to  isolate  the  influence  of  the  tunnel 
boundary  layer  from  the  wing  model  aerodynamics.  The  wing  had  a  chord  and  span  of 
c  =  0.1  m  and  b  =  0.6  m ,  respectively,  producing  a  semi-span  aspect  ratio  of  AR  =  6. 

To  ensure  that  the  structural  dynamics  are  dominated  by  the  virtually  prescribed  tor¬ 
sional  properties,  the  wing  segments  were  constructed  from  molded  urethane  rubber  sections 
with  plastic  ribs  printed  using  a  stereolithography  (SLA)  technique.  The  model  was  built 
in  five  interlocking  sections  allowing  for  quick  modification  of  stiffness  parameters  and  the 
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Figure  2:  Block  diagram  of  the  cyber-physical  control  system.  [Fagley  et  al.,  2015]. 


incorporation  of  different  sensing  transducers  at  various  points  along  the  span.  In  this  con¬ 
figuration  the  model  leading  edge  was  stiffened  by  a  6.4  mm  carbon  fiber  spar  to  ensure  a 
quasi-linear  structure  at  the  tested  wind  speeds,  as  validated  by  the  results  in  Section  4. 
Furthermore,  the  cyber-physical  system  was  used  to  prescribe  a  torsional  stiffness,  torsional 
damping  and  torsional  inertia  of  kg  =  0.47  N  m  rad -1,  r/p  =  6.09  x  10-3  N  m  s  rad ,  and 
Jp  =  4.25  x  10~4  kg  m?  rad _1,  respectively.  These  structural  properties  were  confirmed 
through  system  identification  of  the  model’s  response  to  a  step  input  as  described  in  Fagley 
et  al.  [2015,  2016].  Finally  it  should  be  noted  that  the  stall  flutter  onset  velocity  for  a  nearly 
identical  model  was  experimentally  tested  and  shown  to  occur  at  approximately  20  m/s  in 
Fagley  et  al.  [2016],  placing  the  current  experiments  within  the  post-flutter  regime. 

3.5  Rigid  Wing  Model 

The  test  geometry  for  the  rigid  wing  model  consisted  of  a  finite  span,  NACA  0018  cross- 
section,  rectangular  wing  model  with  a  chord  and  span  of  c  =  0.15  m  and  b  =  0.6  m, 
respectively.  This  produced  a  wing  aspect  ratio,  AR  =  4.  Note  that  the  wing  model  was 
centered  on  a  rotational  shaft  which  symmetrically  spanned  to  the  wind  tunnel  floor  and 
ceiling  producing  a  symmetric  finite  span  wing  with  two  “effectively”  free  tips.  Additionally, 
the  rotation  axis  was  located  at  50%  chord  (i.e.  x/c  =  0.5).  Finally,  the  rigid  wing 
was  constructed  from  three  interlocking  plastic  sections  printed  using  a  stereolithography 
technique. 

For  the  rigid  wing  experiments,  the  cyber-physical  system  dynamics  were  varied  over 
a  range  of  virtually  prescribed  stiffness  and  damping  values  while  the  virtual  inertia,  Jy , 
was  left  fixed  at  zero.  The  physical  interita  Jp  and  physical  mechanical  damping  i]p  of  the 
wing  model  system  were  estimated  through  system  identification  of  the  model’s  response 
to  a  step  input  as  described  in  Fagley  et  al.  [2015,  2016].  This  resulted  in  a  value  of 
Jp  3(10-3)  kg  m2/rad  and  r]p  ps  1.2(10~2)  Nm  s/rad. 

The  rigid  wing  was  also  tested  undergoing  a  prescribed  (or  driven)  sinusoidal  motion  for 
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comparison  with  the  stall-flutter  results  produced  through  the  cyber-physical  control  sys¬ 
tem.  This  was  accomplished  by  switching  the  mechatronic  system  from  the  cyber-physical 
mode  into  a  standard  position  control  mode.  The  position  control  was  accomplished  using 
a  Copely  Controls  Accelnet  Micro  Panel  ACJ-090-12  motor  controller  with  a  0.3  ms  update 
rate  using  a  combination  of  PID  control  with  feed-forward  velocity  and  acceleration  gains 
to  further  reduce  following  errors.  Using  this  approach  the  periodic,  angular  trajectory  of 
the  wing  was  then  defined  by  a  single  frequency  sine  wave  as  given  by  Equation  (2). 

a(t)  =  a  +  Asin(ut)  (2) 

Note  that  the  offset  a ,  the  pitch  amplitude,  A  and  pitch  angular  frequency,  u>,  were  defined 
based  upon  empirical  trajectories  of  the  stall  flutter  limit  cycle  oscillations  measured  in  the 
cyber-physical  mode  of  operation. 

3.6  Stereoscopic  Motion  Capture  System 

The  stereo  vision  data  was  captured  using  two  synchronized  high-speed  Phantom  cameras 
from  Vision  Research  (models  V711  and  V7.3).  The  cameras  were  positioned  in  a  stereo¬ 
scopic  configuration  parallel  to  the  wind  tunnel  floor,  approximately  1.7  m  laterally  away 
from  the  wind  tunnel  center-line  with  a  32°  angular  separation,  as  depicted  in  Figure  3. 
To  provide  the  best  contrast  for  the  motion  capture  program  the  wing  was  painted  black 
and  white  adhesive  dots  were  affixed  to  the  surface.  The  dots  had  a  diameter  of  9  mm  and 
were  positioned  on  the  wing  in  a  rectangular  grid  with  chord-wise  positions  of  x/c  =  0.1, 
0.3,  0.5,  0.7,  and  0.9  and  spanwise  positions  of  z/b  =  0.38,  0.58,  0.78,  and  0.98,  as  detailed 
in  Figure  4.  It  should  be  noted  that  these  spanwise  dot  locations  were  selected  to  align 
with  the  rigid  SLA  ribs  in  the  model.  The  stereo  vision  data  was  captured  at  a  sampling 
frequency  of  fs  =  1  kHz  for  a  period  of  At  =  20  s  during  the  wing  stall  flutter  motion. 

The  image  sequences  were  then  processed  using  a  direct  linear  transformation  (DLT) 
technique  in  the  Matlab  programming  environment  with  the  Digitizing  Tools  version  5 
toolbox  from  Hedrick  [2008] .  The  raw  coordinate  data  produced  by  the  stereo- vision  system 
was  oriented  in  the  standard  wind-tunnel  aligned  Cartesian  system  (x,  y,  z)  as  presented 
in  Figure  5.  Note  that,  the  stereo  vision  measurements  were  only  made  on  the  flexible 
wing  model,  as  a  result  their  discussion  focuses  upon  the  local  twist  displacement,  9,  which 
is  defined  as  the  difference  between  the  instantaneous,  local  angle  of  attack  and  the  root 
angle  of  attack  of  the  flexible  wing  model.  The  translational  displacement  of  the  wing 
coordinates  in  the  cross-stream,  y,  direction  will  also  be  presented  as  Ty.  The  chordwise 
and  spanwise  translational  displacements  of  the  wing  model,  Tx  and  Tz  respectively,  will  not 
be  presented  as  they  were  found  to  be  negligible  and  within  the  experimental  uncertainty. 
The  uncertainties  in  the  cross-stream  bending  and  twist  displacements  were  estimated  to  be 
S(Ty)  =  ±0.85  mm  and  5(6)  =  ±1.6°,  respectively.  Note  that  these  were  determined  from 
a  combination  of  the  calibration  error  and  the  error  associated  camera  pixel  resolution. 
Furthermore,  the  uncertainty  is  computed  from  the  spatial  average  of  the  uncertainties 
across  the  entire  calibration  range  in  the  x  and  y  directions. 

The  stereo-vision  data  was  also  utilized  to  validate  the  tip-twist  measurements  made 
internally  within  the  cyber-physical  system.  As  such  the  stereo-vision  image  capture  was 
also  precisely  synchronized  to  the  measurement  of  the  wing-tip  angular  position  from  the 
motor  shaft  optical  encoder  which  has  a  resolution  of  16,000  counts  per  revolution  or  0.02° 
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Figure  3:  Schematic  defining  the  stereo-vision  camera  orientation  with  respect  to  the  wing 
mounted  in  wind  tunnel  (Top  View). 


z/b  =  0.38  0.58  0.78  0.98 


Figure  4:  Schematic  detailing  the  stereo  vision  tracking  point  layout  along  wing  span  and 
the  coordinate  system  orientation. 


per  count.  For  more  detail  about  the  cyber-physical  system  configuration  see  Fagley  et  al. 
[2015,  2016], 

3.7  Planar  and  Stereoscopic  Particle  Image  Velocimetry 

Both  planar  particle  image  velocimetry  data  (PIV)  and  stereoscopic  particle  image  velocime¬ 
try  (SPIV)  data  was  collected  using  the  same  LaVision  sysem.  The  system  utilized  upto  two 
Imager  SCMOS  cameras  with  2560  x  2160  resolution  and  16  bit  intensity.  The  cameras  were 
fitted  with  either  a  Nikkor  /  =  60  mm  or  a  Nikkor  /  =  105  mm  fixed  focal  length  macro 
lens  producing  scale  factors  of  10.2  pixels/mm  and  19.8  pixels/mm,  respectively.  Note  that 
the  /  =  60  mm  lens  was  used  in  the  planar  PIV  studies  while  the  /  =  105  mm  lenses  were 
used  for  SPIV.  The  flowfield  was  seeded  using  a  ViCount  Compact  1300  smoke  machine 
which  generated  0(1  /jrn)  particles  in  physical  size  or  2  —  3  pixels  in  the  scaled  images. 
The  planar  fields  were  illuminated  by  a  laser  sheet  which  was  produced  by  a  Quantel  Ever¬ 
green  200  532  nrn  dual  pulsed  Nd-YAG  laser  which  had  a  200  mJ  maximum  output  energy. 
For  both  planar  and  stereo  measurements,  the  laser  was  configured  with  a  /  =  —10  mm 
cylindrical  lens  oriented  to  produce  a  laser  sheet  that  illuminated  a  single  chord-wise  sec¬ 
tion  of  the  wing  span.  The  cameras  and  laser  head  were  mounted  on  separate  Velmex  two 
axis  traverses  allowing  for  accurate  positioning  of  the  PIV  system.  The  traverses  had  an 
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Figure  5:  Schematic  defining  the  bending  and  twisting  deflections  in  the  wing  and  wind 
tunnel  coordinate  systems. 


accuracy  of  ±  0.006  mm  and  were  utilized  to  scan  the  wing  Due  to  the  periodic  nature  of 
the  achieved  stall  flutter  oscillation,  the  PIV  data  was  taken  phase  locked  to  the  angular 
position  of  wing  or  flexible  wing  tip,  as  recorded  by  the  cyber-physical  system  encoder. 
Note  that  the  16,  OOOct  encoder  had  a  resolution  of  0.02°. 

For  the  SPIV  study  on  the  flexible  wing  measurements  were  taken  at  multiple  spanwise 
locations,  namely:  z/b  =  0.58,  0.625,  0.71,  0.75,  0.79,  and  0.875.  At  each  spanwise  location, 
three  phase  planes  were  taken  at  <f>  =  204°,  214°,  and  224°  to  capture  the  progression  of 
flow  separation.  Additionally  for  the  spanwise  location  of  z/b  =  0.75,  which  approximately 
corresponded  to  the  location  of  maximum  separation,  higher  resolution  phase  data  was 
collected  at  Acj)  ~  2°  between  <f>  =  194°  and  <fi  =  233°  and  for  two  additional  phase  points  at 
4>  =  185°,  and  243°.  For  each  combination  of  spanwise  and  phase  location,  500  image  sets 
were  collected  with  a  At  =  20  /is.  These  image  sets  were  processed  using  LaVision’s  DaVis 
8.3  program  on  the  graphical  processing  unit  (GPU).  A  multi-pass  scheme  was  utilized 
where  the  initial  pass  utilized  a  guassian  weighted  64  x  64  pixel  interrogation  window  and 
two  subsequent  passes  utilized  32  x  32  pixel  windows,  with  an  overlap  of  75%.  The  raw 
images  were  individually  masked  along  the  surface  reflection  line  resulting  from  the  incident 
laser  light  on  the  wing. 

Planar  particle  image  velocimetry  (PIV)  data  was  also  captured  on  the  rigid  wing  model 
using  the  same  LaVision  system.  In  this  case  a  single  spanwise  location  at  z/b  =  0.5  was 
capture.  Sets  of  25  image-pairs  were  taken  at  64  phase  positions  across  the  full  cycle  for 
the  two  operating  modes.  This  resulted  in  time  steps  of  5  ms  or  phase  steps  of  A</>  =  6°. 
The  planar  PIV  data  was  processed  using  LaVision  DaVis  8.3  software  again  on  the  GPU 
using  a  multi-pass  method  with  an  initial  pass  at  64  x64  pixels  and  two  final  passes  at 
32  x  32  pixels  with  a  50%  overlap  in  the  interrogation  windows. 

3.8  Pitching  Moment  Measurements 

To  understand  the  pitching  moment  results,  a  discussion  of  forces  observed  in  the  cyber¬ 
physical  system  is  required.  The  resulting  moment  response  was  captured  using  a  Futek 
TSS400  single-axis  torque  sensor  installed  inline  between  the  cyber-physical  system  motor 
and  the  wing  on  the  axis  of  rotation.  This  allowed  the  cyber-physical  system  to  measure 
the  net  system  pitching  moment,  TMeas.i  at  the  root  of  the  wing.  This  measured  pitching 
moment  (torque)  represents  the  summation  of  the  aerodynamically  imposed  moment  and 
resistive  moment  imposed  by  the  physical  inertia  and  damping  of  the  wing.  Additionally, 
this  measured  moment  (torque)  is  equal  to  the  torque  produced  by  the  motor.  This  balance 
is  graphically  shown  in  Figure  6  where  r)p  refers  to  physical  damping  imposed  by  the  three 
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Figure  6:  The  balance  of  forces  in  the  rigid-wing  cyber-physical  system. 


bearings  installed  on  the  drive  shaft. 

Mathematically  the  force  balance  is  given  in  Equation  (3): 


T~Motor  —  TM  eas.  —  TAero.  T Struct .  (3) 

where 

T Struct .  —  JpO  +  l]  ['b .  (4) 

The  aerodynamic  pitching  moment,  TAero ■  here  represents  the  resultant  of  all  circula¬ 
tory  and  non-cir dilatory  force  contributions.  This  moment  can  be  re-written  as  a  pitching 
moment  coefficient: 

TAero. 
m  ~  A 

qbc 

which  is  comparable  to  a  traditional  static  moment  coefficient  measure.  In  dynamic  pitching 
situations,  added  mass  appears  as  a  non-circulatory  force  and  depending  on  the  ratio  of  fluid 
density  to  structural  inertia,  this  force  can  be  quite  significant.  In  this  case,  where  the  wing 
operates  in  air,  the  effect  of  this  term  can  be  approximated  by  assuming  the  wing  accelerates 
a  cylinder  of  air  with  radius  c/2  and  length  b.  The  moment  of  inertia  for  this  cylinder  is 
given  by  1/2 Mr2.  Using  this  approximation  the  inertial  resistance  due  to  added  mass  is 
given  by  /07r6c4/32  or  3(1CU5)  kg  m2.  The  wing  system  has  an  inertia  of  ~  3(10~3)  kg  m2 
indicating  that  the  added  mass  is  two  orders  of  magnitude  smaller  and  therefore  negligible. 

The  measured  moment  TMeas  provides  a  net  system  torque  and  is  of  interest  as  it  repre¬ 
sents  the  torque  at  the  wing  root.  This  would  be  the  load  transferred  to  the  aircraft  body,  or 
in  an  energy  harvesting  configuration,  the  net  energy  extracted  from  the  aeroelastic  system. 
This  paper  will  discuss  both  torques  in  terms  of  moment  coefficients,  utilizing  the  following 
two  conventions:  1)  Cm  refers  to  the  moment  coefficient  around  the  50%  chord  location 
resulting  from  the  aerodynamic  pitching  moment  contribution  TAero.  and  2)  CmMeas.  refers 
to  the  moment  coefficient  of  the  measured  or  net  system  torque  TMeas.-  As  a  note,  the  mea¬ 
sured  pitching  moment  TMeas.  is  filtered  using  a  low-pass  butterworth  filter  with  17.5  Hz 
cutoff  frequency. 


4  Results  and  Discussion 

The  results  and  discussion  is  divided  into  two  primary  sections  below.  First,  the  detailed 
investigation  into  the  baseline  kinematics  and  flowfield  development  surrounding  the  flexible 
wing  model  in  stall  flutter  is  presented  in  Section  4.1.  Then,  the  response  comparison  (i.e. 
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Figure  7:  Time-history  of  the  local  angle  of  attack,  a,  obtained  by  the  cyber-physical 
tip-twist  encoder  (black)  compared  against  the  stereo-vision  motion  tracking  system  mea¬ 
surements  at  spanwise  positions  of  z/b  =  0.38,  0.58,  0.78,  and  0.98  (magenta,  yellow,  red, 
and  blue  respectively). 


kinematic  and  flowfield)  is  made  in  Section  4.2  for  the  rigid  wing  model  undergoing  both 
passive  stall-flutter  limit  cycle  oscillations  and  prescribed  pitch  oscillations. 

4.1  Flexible  Wing  Model 

The  resulting  structural  kinematics  of  the  high-aspect-ratio,  finite-span,  flexible,  cyber¬ 
physical  wing  system  undergoing  stall  flutter  are  thoroughly  presented  in  five  subsections 
sections  below.  First  in  Section  4.1.1,  a  comparison  is  made  between  the  cyber-physical 
system’s  encoder  tip-twist  angle  measurement  and  the  twist  measurements  made  via  the 
stereo- vision  system;  to  cross- validate  both  measurement  systems.  Second  in  Section  4.1.2, 
a  detailed  discussion  is  presented  of  the  trajectory,  periodicity,  and  uniformity  of  the  stall 
flutter  limit  cycle  oscillations.  Third,  a  fourier  analysis  of  the  time  history  is  conducted 
in  Section  4.1.3  which  then  leads  to  the  eigenfunction  projection  in  Section  4.1.4.  The 
time-signals  for  each  of  the  principal  motions,  torsion  and  bending,  are  decomposed  into 
the  modes  given  by  the  eigenfunction  solutions  from  classical  beam  vibration  theory.  From 
this  projection  and  fourier  analysis,  a  kinematic  relation  for  the  wing  motion  is  written  in 
Section  4.1.5.  Finally,  following  the  discussion  of  the  wing  kinematics,  the  aerodynamics 
of  the  flexible  wing  system  are  thoroughly  discussed.  This  begins  with  the  analysis  of  the 
internal  moments  (torques),  measured  from  the  cyber-physical  system,  in  Section  4.1.6. 
Then  the  flowfield  surrounding  the  flexible  wing  is  analyzed  and  discussed  in  Sections  4.1.7, 
4.1.8,  and  4.1.9  where  the  both  the  periodic  and  spatial  variations  in  the  flowfield  discussed. 

4.1.1  Validation  and  Characterization  of  Stereo-Vision  Measurements 

In  order  to  ensure  that  the  stereo- vision  system  is  performing  as  expected,  a  simple  cross- 
validation  with  the  tip  encoder  is  conducted.  Figure  7  presents  the  time- history  of  the  local 
angle  of  attack  at  four  spanswise  locations  measured  by  the  stereo-vision  system  plotted 
with  the  wing  tip  angle  of  attack  measured  by  the  cyber-physical  system’s  optical  encoder. 
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From  these  results  a  very  close  agreement  is  observed  between  the  cyber-physical  predicted 
tip  angle  of  attack  (tip-twist  angle)  and  the  stereo-vision  angle  of  attack  at  98%  span.  The 
signals  match  exceptionally  well  in  phase,  but  display  small  deviations  in  the  peak  amplitude 
up  to  Aa  =  2°.  While  this  deviation  is  on  the  order  of  the  measurement  uncertainty  in  the 
stereo-vision  system  as  reported  in  Section  3.6,  it  appears  to  be  a  systematic  error  in  the 
experiment.  Specifically,  the  stereo-vision  data  consistently  under- measures  the  peak  angle 
of  attack  in  the  positive  direction,  while  consistently  over-measuring  the  negative  peak. 
This  suggests  that  there  is  a  small  bias  between  the  a  =  0°  positions  in  the  cyber-physical 
and  stereo-vision  measurement  systems  on  the  order  of  (5(cc)  ~  0.5°.  Additionally  the  larger 
discrepancy  at  the  positive  peaks  may  also  be  due  to  the  exacerbation  of  the  stereo-vision 
measurement  uncertainty  with  greater  distance  away  from  the  cameras.  Specifically,  note 
that  there  is  similar  phasing  between  the  twisting  and  bending  modes  of  the  wing;  thus 
as  the  wing  reaches  its  maximum  angle  of  attack,  axip  =  30°,  the  surface  also  reaches  its 
maximum  bending  deflection.  All  that  said,  the  agreement  validates  the  captured  stereo¬ 
vision  data  within  the  confidence  of  the  measurement  systems. 

At  the  given  testing  conditions,  noted  in  Section  3.4,  the  cyber-physical  wing  system 
underwent  a  very  periodic  oscillation.  This  is  easily  observed  through  examination  of  the 
cyber-physical  data  and  the  processed  stereo-vision  motion  as  was  shown  in  Figure  7.  The 
tip  notably  followed  a  sinusoidal-like  motion  with  a  frequency  of  /  =  5.4  Hz  and  an  ampli¬ 
tude  of  approximately  15°  centered  at  axip  =  14.5°.  It  should  be  noted  that  the  wing  root 
angle  of  attack  was  fixed  at  aR00t  =  10°,  thus  the  tip  on  average  displayed  a  5°  washout 
relative  to  the  root  position. 

Further  visual  examination  of  this  limit  cycle  oscillation  presented  in  Figure  7  shows  that 
the  pitching  oscillation  is  dependent  on  spanwise  position.  This  is  indicative  of  the  spanwise 
gradients  expected  from  the  flexible  (or  deformable)  wing  model.  It  is  also  apparent  that 
the  pitching  motion  is  biased  towards  positive  angles  of  attack.  Specifically,  the  pitch  angle 
oscillates  from  8  =  —5°  to  8  ~  30°  with  respect  to  the  root  angle  of  attack  (i.e.  the  root 
referenced  twist  angle). 

4.1.2  Periodicity  and  Cycle  Variation  in  the  Wing  Deformation 

Before  discussing  the  modal  analysis,  the  variability  in  the  stall  flutter  oscillation  is  under¬ 
stood.  To  visualize  this  variability,  a  cycle  to  cycle  overlay  of  the  pitching  and  the  bending 
signals  at  the  tip  location  across  the  entire  20  s  period  is  provided  in  Figure  8.  Addition¬ 
ally  in  Figure  8a,  a  sine  wave  is  plotted  for  reference  on  top  of  the  tip-twist  data  (dashed 
red  line).  This  clearly  demonstrates  a  noticeable  distortion  of  the  trajectories  from  a  pure 
sinusoidal  motion.  For  the  tip-twist  data,  this  results  qualitatively  in  the  sharpening  and 
broadening  of  the  positive  and  negative  peaks  in  the  cycle  respectively.  The  influence  of 
the  secondary  frequency  is  even  more  apparent  in  the  bending  overlay,  where  secondary 
oscillations  appear  to  be  superposed  on-top  of  the  primary  frequency.  From  this  overlay, 
the  variability  in  amplitude  is  qualitatively  observed  to  be  on  the  order  of  A 8  ~  ±3°  and 
ATy  <  ±10  mm. 

The  cycle  to  cycle  variability  in  frequency  or  phase  is  demonstrated  when  taking  the 
autocorrelation  of  the  raw  x ,  y  displacement  data.  Figure  9  presents  the  normalized,  au¬ 
tocorrelation  in  time  for  each  of  the  torsional  and  bending  displacement  time-series,  and 
strongly  demonstrates  the  wandering  of  the  kinematics  through  the  testing  window.  Each 
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Figure  8:  Cycle-to-cycle  variability  present  in  data  as  shown  by  overlaying  each  individ¬ 
ual  cycle  of  the  20  s  data  capture  period  (approximately  105  cycles)  for  (a)  the  torsional 
deflection,  9,  with  sine  wave  plotted  for  reference,  and  (b)  the  cross-stream  bending,  Ty. 


of  the  displacements  remain  highly  correlated  within  the  first  two  seconds  of  the  full  20  s 
period,  implying  that  the  LCO  is  stable  and  highly  periodic  for  approximately  11  cycles. 
However  after  the  initial  2  s  a  significant  reduction  in  correlation  value  is  observed,  indi¬ 
cating  that  the  kinematics  are  wandering  away  from  the  initial  state.  At  approximately  the 
4  s  mark  the  correlation  drops  to  near  zero  values.  Since  the  cycle-to-cycle  overlay  data 
presented  previously  in  Figure  8,  does  not  show  a  very  large  change  in  waveform  shape, 
it  is  inferred  that  this  decrease  in  correlation  represent  changes  in  frequency  or  phase  and 
amplitude  as  time  progresses. 

Figure  10  shows  phase  plane  plots,  where  the  bending  deflection  is  plotted  against  the 
torsional  deflection  for  three  limit  cycles.  This  data  shows  that  the  wing  motion  clearly 
follows  a  repeated,  limit  cycle  trajectory  with  a  significant  hysteresis  that  grows  uniformly 
with  spanwise  position.  It  also  demonstrates  the  wandering  in  the  limit  cycle  trajectory 
that  was  previously  discussed.  As  a  note,  plotting  over  more  cycles  would  also  show  an 
increase  in  variability  for  this  figure. 

From  this  data,  the  maximum  phase- averaged  cross-stream  bending  deflection,  Ty,  ap¬ 
pears  to  correlate  closely  with  the  maximum  pitch  angle,  6 ,  at  each  spanwise  position. 
Running  a  cross-correlation  indicates  a  phase  lead  in  the  bending  signal  of  A t/T  =  0.032 
(or  a  phase  angle  of  7  =  11.5°  in  the  average  cycle).  This  can  potentially  be  explained 
when  considering  the  aerodynamic  force  production.  The  largest  lift  force  will  occur  near 
the  maximum  angle  of  attack  just  before  the  flow  field  fully  separates  and  the  wing  stalls. 
Additionally,  the  wing’s  bending  response  time  is  expected  to  be  much  faster  than  its  tor¬ 
sional  response  time  owing  to  its  lower  torsional  natural  frequency  as  compared  with  bend¬ 
ing.  Thus  it  is  inferred  that  the  maximum  bending  deflection  correlates  closely  to  the  stall 
point,  but  that  the  torsional  motion  persists  beyond  the  stall  point.  This  is  potentially  due 
to  the  lower  torsional  stiffness  of  the  system,  but  could  also  be  related  to  cycle  differences  in 
the  maximum  lift  and  pitch  moment  production;  which  requires  further  investigation.  Even 
so  the  structural  LCOs  are  primarily  dominated  by  the  stall  aerodynamics  of  the  wing;  as 
previously  observed  with  stall  flutter  [Dimitriadis  and  Li,  2009,  Razak  et  al.,  2011,  Tang 
and  Dowell,  2001,  2002], 

Finally,  it  should  be  noted  that  the  wing  was  designed  to  be  torsionally  weak,  but  strong 
in  bending  to  effectively  decouple  the  modes.  As  shown  in  the  preceding  results,  the  bending 
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(a)  Torsion 


t,  s 

(b)  Cross-stream  Bending 

Figure  9:  Autocorrelation  in  time  of  the  (a)  torsional  displacement,  6  and  (b)  cross-stream 
bending  displacement,  Ty  signals. 


signal  is  not  fully  decoupled.  This  is  because  the  model  has  finite  stiffness  in  bending  and 
as  discussed  previously,  the  aerodynamics  produce  a  maximum  lift  force  near  the  maximum 
pitch  angle.  With  that  said,  the  2  cm  bending  deflection  is  quite  small  in  comparison  to 
the  30°  pitching  amplitude,  indicating  that  the  wing  is  dominated  by  the  torsional  motion, 
as  intended. 

4.1.3  Fourier  Decomposition  of  the  Wing  Deformation 

Evidence  of  a  single  (or  first)  mode  dominance  appears  in  the  spectral  analysis  of  the 
pitching  and  bending  signals  presented  in  Figures  11  and  12  where  a  fast  Fourier  trans¬ 
form  (FFT)  was  applied  to  the  torsion  and  cross-stream  bending  for  each  of  the  spanwise 
positions  independently.  Note  that  the  number  of  samples  used  in  the  FFT  was  16384 
(T  ~  16  -s)  at  a  sampling  frequency,  fs  =  1  kHz  thus  the  frequency  resolution  in  the  plots 
are  A /  =  0.061  Hz.  In  both  of  the  plots,  a  dominant  frequency  of  /  =  5.4  Hz  and  sec¬ 
ond  frequency  at  /  =  10.7  Hz  (approximately  twice  that  frequency)  are  observed  in  both 
signals  though  there  is  some  broadening  of  the  FFT  peaks.  This  is  likely  a  result  of  the 
frequency  variation  discussed  previously.  The  dominant  frequency  matches  closely  to  the 
expected  torsional  natural  frequency  of  the  flexible  wing  system,  which  is  estimated  to  be  to 
be  ojQ  =  5.66  Hz  using  the  prescribed  properties  of  the  cyber-physical  system  (Section  4.1) 
and  Equation  7  for  a  second-order  linear  dynamic  system.  The  observed  second  frequency, 
/  =  10.7  Hz ,  distorts  the  torsional  motion  away  from  a  pure  sinusoidal  oscillation,  inducing 
the  peak  sharpening/broadening  observed  in  Figure  8a.  Furthermore,  no  oscillations  are 
observed  around  the  primary  frequency  since  the  secondary  frequency  is  only  twice  that  of 
the  primary. 
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Figure  10:  Phase-plane  plot  of  Ty  bending  deflection  versus  torsional  deflection  angle,  9 ,  at 
spanwise  locations  of  z/b  =  1,0.78,0.58,  and  0.38,  over  three  limit  cycle  periods. 
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Plotted  alongside  the  frequency-amplitude  spectrum  for  each  signal  at  each  spanwise 
location  is  a  frequency-phase  plot.  In  the  torsion  signal,  this  plot  shows  that  the  first 
frequency  components  have  matching  phase  angles  while  the  second  frequency  component 
shows  a  significant  variation  in  phase.  The  phase  locking  of  the  dominant  frequency  indicates 
that  the  wing  motion  for  that  frequency  can  be  described  by  a  single,  continuous  mode  shape 
oscillating  with  the  same  sinusoidal  trajectory.  The  variability  in  the  calculated  phase 
angle  for  the  second  frequency  is  more  likely  an  indication  of  data  noise  as  the  associated 
amplitudes  for  this  fourier  component  are  quite  small,  though  as  will  be  shown  later  they  are 
still  significant.  Subsequently,  all  wing  sections  are  considered  to  move  torsionally  in-phase. 

Examining  the  phase  of  the  bending  signal  indicates  two  frequencies  are  also  present 
however  they  share  similar  dominance.  Furthermore,  the  phase  angles  for  both  of  the 
frequencies  match  closely;  which  may  further  enforce  the  argument  that  the  phase  variation 
in  the  second  torsional  frequency  is  more  a  result  of  small  amplitude  signals  than  a  physical 
effect.  Also  note  in  all  of  these  decompositions,  the  fourier  spectrum  shows  a  significant 
non-zero  mean  which  validates  that  the  wing  oscillates  around  a  position  different  from 
the  base  angle-of-attack.  Additionally,  since  this  mean  offset  is  positive,  the  wing  favors 
positive  deflection. 

Finally,  it  should  also  be  noted  that  the  bending  natural  frequency  of  a  nearly  identical 
flexible  wing  system  was  measured  to  be  u>h  =  11.45  Hz  by  Fagley  et  al.  [2016].  This 
frequency  closely  matches  the  second  dominate  frequency  peak  observed  in  the  data  and 
may  also  explain  why  this  frequency  appears  stronger  in  the  bending  spectrum.  That 
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Figure  11:  Comparison  of  a)  amplitude  and  b)  phase  of  spectral  (FFTs)  analysis  of  Torsion 
signal  projected  to  rigid  body. 
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Figure  12:  Comparison  of  a)  amplitude  and  b)  phase  of  spectral  (FFTs)  analysis  of  cross 
stream  bending  signal  projected  to  rigid  body. 
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said  the  first  torsional  frequency  still  dominates  the  spectral  content,  with  no  observable 
frequencies  related  to  the  higher  order  torsion  or  bending  eigenmodes.  Note  that  these  would 
theoretically  be  expected  to  appear  at  wemode2  =  3a ig  and  Uhmode2  =  6.3a Jh,  respectively  for  a 
cantilevered  beam  in  free  vibration  [Hodges  and  Pierce,  2014].  It  is  however  possible,  that 
due  to  spatial  aliasing  across  the  wing,  these  frequencies  are  not  observed  in  the  FFT. 

4.1.4  Eigenfunction  Projection  of  the  Wing  Motion 

From  the  Fourier  expansion  it  is  expected  that  a  single  dominant  mode  (in  the  frequency 
domain)  will  appear  in  the  torsional  motion,  potentially  oscillating  sinusoidally.  In  contrast 
for  the  bending  motion  two  equally  dominant  modes  are  expected.  Therefore  to  provide  a 
more  descriptive  decomposition  of  the  wing  motion,  in  the  time  domain,  the  eigenfunction 
solutions  for  a  cantilevered  beam  are  utilized.  This  was  accomplished  through  projecting 
the  measured  unsteady  torsion  and  bending  deflections  onto  the  first  three  analytically 
predicted  mode  shapes  for  a  clamped-free  beam  in  both  torsion  and  bending  vibration. 
While  the  Fourier  expansion  doesn’t  show  frequency  peaks  associated  with  higher  modes, 
three  eigenmode  shapes  are  utilized  to  ensure  that  the  non-linearities  in  the  wing  motion 
can  be  best  captured.  The  unsteady  (time-varying)  mode  weights  associated  with  these 
modal  projections  are  plotted  in  Figures  13  and  14  for  the  torsional  component,  6,  and  the 
cross-stream  bending  component,  Ty.  respectively.  Note  that  the  Fourier  decompositions  of 
each  of  the  unsteady  mode  weights  are  also  plotted  in  these  figures. 

The  torsional  time-varying  mode  weights,  Figure  13a,  similarly  display  the  dominance 
of  the  first  torsional  mode  in  the  wing  motion.  Specifically,  the  first  mode  displays  the 
largest  amplitude,  as  compared  with  the  amplitudes  of  the  second  and  third  modes  which 
are  nearly  an  order  of  magnitude  smaller.  Therefore,  the  majority  of  the  motion  is  derived 
from  the  first  eigenmode  with  the  second  and  third  potentially  capturing  the  small,  time 
varying  perturbations  around  this  principle  mode. 

As  with  the  torsion  modes,  Figure  14a  presents  the  time- varying  bending  mode  weights. 
The  bending  deflection  shows  a  similar  dominance  in  the  first  mode,  however  the  second 
and  third  modes  show  more  significant  weight  than  was  observed  in  the  torsional  motion. 
The  bending  deflection  also  shows  a  changing  amplitude  through  the  20  s  window  that  may 
be  consistent  with  a  low  frequency  beat.  It  is  reasonable  to  expect  a  beat  to  form  in  the 
bending  signal,  because  as  noted  previously,  the  second  frequency  in  the  FFT  spectrum  is 
very  near  to  the  expected  natural  frequency  of  the  wing  in  bending. 

As  expected,  the  spectral  content  of  the  modes,  in  each  the  torsion  and  bending  direc¬ 
tions,  matches  closely  to  that  of  the  data  prior  to  projection  onto  the  theoretical  eigenmodes 
(Figures  11  and  12).  Furthermore,  the  first  mode  clearly  dominates  the  frequency  content, 
though  the  second  mode  does  have  greater  importance  in  the  bending  signal. 

This  analysis  demonstrates  that  within  the  20  s  windows,  the  waveform  shape  is  ac¬ 
curately  represented  by  the  eigenmode  decomposition  which  is  also  able  to  capture  the 
wandering  of  the  signals  over  time.  Since  the  loss  in  correlation  is  attributed  to  small  fre¬ 
quency  and  amplitude  variations,  as  discussed  in  Section  4.1.2,  rather  than  fundamental 
changes  in  waveform  shape,  these  modes  can  be  generalized  and  expected  to  describe  the 
motion  for  an  arbitrary  window  size.  Most  notably,  the  linear  mode  shapes  accurately 
capture  the  nonlinear  wing  response  in  both  bending  and  torsion. 

Taking  this  approach,  the  non-linearities  appear  in  the  time-varying  modal  weights 
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(FFTs)  analysis  of  Torsion  eigenmode  signal  for  first  three  modes. 
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Figure  14:  Comparison  of  a)  bending  mode  weights  b)  amplitude  and  c)  phase  of  spectral 
(FFTs)  analysis  of  bending  eigenmode  signal  for  first  three  modes. 
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associated  with  each  mode.  However,  with  the  dominance  of  the  first  mode,  in  the  torsional 
deflection  and  the  relatively  small  bending  deflection,  the  wing  deformation  can  be  largely 
constructed  by  the  first  mode  in  torsion  and  its  associated  weight.  For  completeness,  the 
next  section  will  determine  weights  associated  with  each  of  the  modal  components. 

4.1.5  Kinematic  Equations 

Based  upon  the  eigenmode  projection  and  the  fundamental  frequencies  present  in  the  data, 
an  empirically  defined  kinematic  relation  can  be  written  for  the  stall  flutter  motion  observed. 

In  the  torsional  direction  (Figure  13),  the  second  frequency  at  10.7  Hz  is  relatively 
small  in  amplitude  as  compared  to  the  first  (dominate)  frequency,  however  its  presence  is 
still  significant.  Specifically,  it  was  previously  shown  that  the  motion  is  not  purely  sinusoidal 
and  thus  the  second  frequency  influences  the  asymmetry,  or  peak  broadening/sharpening, 
across  the  cycle.  Additionally, the  higher  order  modes  are  retained  to  potentially  capture 
spanwise  non-linearities  in  the  wing.  Also  note  that  inspection  of  the  phase  plot  (Figure  13c) 
shows  that  the  second  and  third  modes  are  nearly  180°  out  of  phase  for  the  torsional 
component.  However  this  does  not  imply  cancellation  of  these  modes,  as  they  represent 
orthogonal  shapes  whose  spatial  dependencies  are  also  180°  out  of  phase. 

In  contrast,  the  Fourier  decomposition  of  the  bending  mode  weights  (Figure  14)  shows 
that  the  two  frequencies  share  similar  levels  of  dominance  in  both  the  first  and  second  modes 
of  bending.  However  the  third  mode  shows  no  significant  spectral  content,  thus  it  will  be 
ignored  in  the  mathematical  model  of  the  bending  motion. 

The  kinematic  equations  of  motion  for  the  wing  are  thus  written  in  terms  of  both  the 
eigenmodes,  which  prescribe  the  spatial  dependencies,  and  the  Fourier  series  components  at 
the  identified  frequencies,  which  prescribe  the  time  dependence.  The  resulting  relationships 
are  defined  by  Equations  8  and  9.  These  equations  are  consistent  with  previous  work  in 
that  they  isolate  the  motion  to  a  principally  torsional  oscillation  (one  degree  of  freedom), 
but  allow  for  small  bending  components  to  couple  in.  In  contrast  to  previous  works,  the 
kinematic  relation  shows  the  presence  of  two  frequencies  indicating  that  the  stall  flutter 
motion  is  not  restricted  to  a  pure  sine  wave  in  a  deformable  wing  body. 

2 

+  y ^[aijcos 

3= 1 


For  the  current  experimental  data,  all  of  the  empirical  coefficients  were  determined 
through  a  least  squares  fit  across  the  full  20-s  period  of  the  data  set.  Note  that  the  fun¬ 
damental  frequencies  were  found  to  be  uq  =  5.4  Hz  ±0.20  Hz  (torsion)  and  uq  =  5.4  Hz 
±0.59  Hz  (bending)  and  uq  =  10.7  Hz  ±0.31  Hz  (torsion)  and  uq  =  10.7  Hz  ±0.59  Hz 
(bending).  The  remaining  coefficients  are  included  in  Tables  1  and  2  for  the  torsion  and 
bending  components,  respectively.  The  mean  coefficients  as  well  as  the  standard  deviations 
from  these  values  are  included  in  the  tables.  Also  note  that  the  phase  shift  between  the 
bending  and  torsional  signals  is  captured  in  the  variable  7,  and  was  found  to  be  7  =  11.5°, 
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Table  1:  Empirical  coefficients  for  the  torsional  kinematics,  0(z,t)  (Equation  8). 
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Table  2:  Empirical  coefficients  for  the  cross-stream  translation  kinematics,  Ty(z,t)  (Equa¬ 
tion  9). 


WTy,  i 

1.2  ±0.65 

Wtv,2 

0.26  ±0.054 

Cll 

-0.57  ±1.7 

C21 

-0.09  ±0.10 

Cl  2 

-1.0  ±0.76 

C22 

-0.19  ±0.09 

dn 

-2.3  ±1.8 

02  i 

-0.16  ±0.12 

d\2 

0.33  ±1.06 

022 

0.023  ±0.10 

with  the  torsional  signal  lagging  the  bending.  The  tabulated  values  in  Tables  1  and  2  are 
consistent  with  previous  arguments  that  the  pitching  motion  is  dominated  by  a  single  tor¬ 
sional  mode  with  slight  perturbations  from  higher  order  modes  as  demonstrated  by  their 
small  magnitudes. 

The  mean  wing  trajectory  is  reconstructed  and  compared  to  the  raw  cycle  to  cycle  data 
for  the  full  20s  period  in  Figures  15  and  16  utilizing  the  empirically  determined  kinematic 
relations  for  the  torsional  (Equation  8)  and  bending  (Equation  9)  motions,  respectively. 
Note,  that  for  the  torsional  component  (Figure  15)  a  progression  of  modal  reconstructions 
are  presented  including:  (1)  only  the  dominate  (#i),  (2)  the  first  two  (0 1  +  0-2 ),  and  (3)  all 
three  eigenmodes  ( 01+02+03 )•  All  three  reconstructions  of  the  torsional  motion  qualitatively 
capture  the  cycle  to  cycle  mean  to  approximately  80%  for  the  full  model  span.  That 
said,  it  is  clear  that  utilizing  only  the  first,  dominate  eigenmode  causes  the  tip  amplitude 
( z/b  =  98%)  to  be  under-predicted,  while  significantly  over-predicting  the  amplitude  at 
z/b  =  38%.  Progressively  including  the  additional  eigenmodes  allows  the  model  to  better 
capture  the  spanwise  non-linearity,  but  with  diminishing  returns. 

Greater  deviations  from  the  raw  data  are  visible  in  the  reconstruction  of  the  mean 
bending  signal  (Figure  16).  This  may  be  the  result  of  a  greater  variation  between  the 
cycles  across  the  data  period,  which  would  skew  the  mean  coefficients.  With  that  said, 
the  coefficients  can  be  individually  fit  to  any  individual  cycle  allowing  for  a  nearly  exact 
recreation  of  that  cycle.  However  further  work  would  be  required  to  better  adapt  this 
approach  to  a  fully  deterministic,  predictive  model  for  the  entire  time  history  (20  s  period). 
More  specifically,  a  model  that  accurately  describes  the  time  dependent  phase  and  amplitude 
variability  across  the  data  set  would  be  required.  Thus  the  primary  limitation  of  the  current 
approach  lies  within  the  independent  (non-linear)  weightings  of  each  eigenmode.  While  this 
allows  the  empirical  model  to  be  adapted  to  non-linear  variations  in  the  wing  motion,  it 
limits  the  predictive  ability  of  the  model.  With  that  said  the  results  show  that  the  torsional 
deflection  can  be  predicted  across  the  span  to  approximately  80%  accuracy  by  the  first 
torsional  mode  alone.  Thus  for  some  in-situ  applications  it  may  be  sufficient  to  predict  the 
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Figure  15:  Comparison  of  the  raw  and  reconstructed  trajectories  in  6  for  each  spanwise 
measurement  location. 


spanwise  deformation  from  only  the  first  mode. 

4.1.6  Integrated  Torque  Prescribed  by  Cyber-Physical  System 

To  gain  some  insight  into  the  development  of  forces  around  the  flexible,  finite  span  wing, 
the  moment  developed  at  the  wing  tip  due  to  the  flowfield  was  calculated  using  the  method 
discussed  in  Section  3.8.  As  was  discussed  previously,  this  tip  torque  is  not  the  aerodynamic 
pitching  moment  as  would  be  measured  at  the  wing  root,  but  is  rather  related  to  the 
net  deflection  of  the  wing  tip  resulting  from  all  aerodynamic  pitching  moments  integrated 
along  the  surface  of  the  wing.  This  tip  moment  is  normalized  and  presented  as  a  moment 
coefficient  using  Cm  =  Tup/qSc.  In  the  proceeding  analysis,  two  data  cases  are  analyzed 
where  the  first  corresponds  to  the  stereo- vision  data  published  in  Culler  et  al.  [2017],  where 
kg  =  0.47  N  m  rad -1,  r]  =  6.09  x  10-3  N  m  s  rad -1,  and  J  =  4.25  x  10-4  kg  m2  rad -1. 
The  second  corresponds  to  the  structural  properties  during  the  SPIV  testing  and  are  kg  = 
0.54  N  m  rad -1,  rj  =  5.0  x  10-3  N  m  s  rad-1,  and  J  =  4.25  x  10-4  kg  m 2  rad-1.  These 
cases  will  be  subsequently  referred  to  as  the  stereo-vision  and  SPIV  cases,  respectively.  In 
both  data  sets,  the  same  wing  model  was  tested  at  the  same  tunnel  velocity  of  U0 0  =  26m/s. 
Because  the  same  model  was  utilized  and  the  virtual  inertia  Jy  was  left  constant,  the  inertia 
was  the  same  for  each  test  case.  The  stiffness  and  damping  were  varied  slightly  between 
the  cases  by  0.07  N  m  rad-1  and  1.09  x  10-3  N  m  s  rad-1  respectively,  where  the  first 
case  had  lower  stiffness  but  higher  damping  and  the  second  had  higher  stiffness  but  lower 
damping.  These  structural  properties  were  estimated  through  system  identification  of  the 
model’s  response  to  a  step  input  as  described  in  Fagley  et  al.  [2015]. 

The  calculated  Cm  is  plotted  against  the  tip  trajectory  in  Figure  17  for  both  cases, 
showing  a  clear  hysteresis  loop.  The  figure  shows  that  both  Cm  versus  6  loops  have  similar 
character  however  the  case  corresponding  to  the  SPIV  data  (red)  reached  a  negative  angle 
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Figure  16:  Comparison  of  the  raw  and  reconstructed  trajectories  in  Ty  for  each  spanwise 
measurement  location. 


of  attack  approximately  4°  lower  than  the  stereo-vision  case  (grey).  This  is  consistent 
with  the  lower  structural  damping  in  the  SPIV  case  which  is  expected  to  allow  for  larger 
amplitude  pitch  oscillations.  In  both  cases,  the  pitch-up  portion  of  motion  was  characterized 
by  positive  torques  and  the  pitch  down  motion  by  small,  negative  net  torque.  This  implied 
that  the  integrated  aerodynamic  pitching  moment  along  the  wing  largely  goes  to  zero  and 
that  the  wing  angular  position  is  restored  primarily  due  to  the  internal  structural  forces. 
This  further  suggests  that  the  wing  was  fully  stalled  during  the  pitch  down  portion  of  the 
motion.  This  significant  decay  in  the  pitching  moment  is  also  consistent  with  the  general 
expectations  of  dynamic  stall,  where  the  lift  force  significantly  drops  after  a  stall  condition 
is  reached.  It  should  also  be  noted  that  the  wing  reached  a  maximum  tip  angle  of  attack 
of  OLmax  ~  30°  ( 6max  ~  20°)  for  both  the  stereo-vision  and  SPIV  cases.  An  angle  of  attack 
of  30°  is  significantly  higher  than  the  static  stall  angle  of  19°  as  presented  in  Goett  and 
Bullivant  [1938]  and  reaching  this  high  angle  of  attack  before  stall  is  consistent  with  classic 
descriptions  of  dynamic  stall  [Larsen  et  al.,  2007,  Leishman  and  Beddoes,  1989]. 

From  the  calculation  of  internal  pitching  moment,  the  stability  of  the  coupled  system 
can  be  determined  through  calculation  of  differential  aerodynamic  work  as  discussed  byM- 
cCroskey  [1981].  Specifically,  the  direction  of  the  differential  work,  dW  =  — Cmda  will 
indicate  whether  the  structure  is  doing  work  on  the  flow  or  whether  the  flow  is  doing  work 
on  the  structure.  In  the  case  of  positive  dW  the  structure  is  doing  work  on  the  flow  and  the 
system  is  stable.  In  the  case  of  the  negative  dW  the  fluid  is  doing  work  on  the  structure 
and  the  system  will  tend  towards  instability,  or  growth  of  LCO  amplitude. 

A  plot  of  dW  versus  tip-twist  angle,  9,  is  provided  in  Figure  18.  The  figure  shows  that 
the  peak  negative  dW  is  approximately  17%  larger  for  the  SPIV  case  (red)  than  for  the 
stereo- vision  case  (grey).  This  is  again  consistent  with  the  lower  structural  damping,  which 
creates  a  less  stable  configuration.  From  the  figure,  it  is  further  observed  that  the  dW 
becomes  positive  at  the  turning  points  or  maximum  and  minimum  twist  angles.  This  also 
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Figure  17:  Cm  versus  tip-twist  angle,  9 ,  loops  for  the  stereo- vision  case  (grey)  and  the  SPIV 
case  (red). 
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Figure  18:  The  differential  work  dW  versus  tip-twist  angle,  9  for  the  stereo-vision  case 
(grey)  and  the  SPIV  case  (red). 


makes  sense  as  it  indicates  that  the  stabilization  of  the  system  provides  the  bounds  on  the 


amplitude  of  oscillation.  From  the  dynamic  stall  perspective  as  discussed  previously,  this 
maximum  amplitude  boundary  is  likely  associated  with  wing  stall.  Furthermore,  during  the 
majority  of  the  cycle,  dW  is  negative  indicating  that  the  aerodynamics  are  dominating  the 
structural  response  as  was  observed  by  Tang  and  Dowell  [2001,  2002], 


From  these  results  it  becomes  clear  that  the  observed  stall  flutter  LCO  exhibits  a  delay 


of  stall  behavior  that  is  likely  attributable  to  a  dynamic  stall  event.  Furthermore,  the  results 
show  that  the  pitching  amplitude  maxima  is  associated  with  a  re-stabilization  of  the  system. 
This  stability  bound  is  likely  associated  with  stall  of  the  wing.  With  this  information,  stereo 
particle  image  velocimetry  data  was  collected  near  the  peak  angle  of  attack,  encompassing 
the  phase  region  where  flow  separation  occurs.  These  phases  where  selected  to  determine 
whether  a  dynamic  stall  vortex  develops  and  interacts  with  the  wing  surface  as  is  the  case 
in  rigid  body  motion. 
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(a)  Angle  of  Attack  (b)  Pitching  Moment 

Figure  19:  Phase  locations  of  SPIV  data  collection  relative  to  the  (a)  wing  tip  angle  of 
attack  and  (b)  internal  pitching  moment. 

4.1.7  Flow  Field  Measurements 

In  order  to  understand  how  flow  separation  evolves  along  the  wing,  stereoscopic  particle 
image  velocimetry  measurements  where  taken  at  the  spanwise  and  phase  locations  noted 
in  Section  3.7.  Specifically,  a  phase  sweep  at  z/b  =  0.75  is  presented,  corresponding  to  the 
approximate  spanwise  location  of  maximum  separation.  SPIV  planes  at  additional  spanwise 
locations  of  z/b  =  0.58,  0.625,  0.71,  0.75,  0.79,  and  0.875  are  presented  to  determine  whether 
or  not  the  stall  flutter  flowfield  was  characterized  by  widespread  or  localized  separation. 
All  of  the  captured  phase  locations  for  each  spanwise  position  are  located  within  the  region 
illuminated  in  Figure  19. 

As  the  figure  shows,  the  SPIV  phases  start  at  the  approximate  peak  in  Cm  and  progress 
through  the  region  of  decreasing  Cm.  As  a  note,  despite  the  general  trend  of  decreasing  Cm, 
all  of  the  SPIV  phases  correspond  to  pitch  up  motion  of  the  wing,  as  also  demonstrated 
in  Figure  19.  Thus  the  maximum  Cm  occurs  well  before  the  wing  reaches  maximum  pitch 
angle. 

The  results  for  the  phase  sweep  at  z/b  =  0.75  are  presented  in  Section  4.1.8  while 
the  results  for  the  spanwise  sweep  are  presented  in  Section  4.1.9.  For  the  phase  sweep  at 
z/b  =  0.75  select  phases  are  presented,  interspersed  across  the  measurement  window,  to 
show  the  principle  features  of  separation  for  the  flexible  wing. 

4.1.8  Phase  Evolution  of  the  Flow  Field 

Figure  20  presents  streamlines  overlaid  on  normalized  vorticity  contours,  ilz  jj—  for  phases 
angles  of  cj)  =  185°,  200°,  216°,  224°,  231°,  and  243°  at  z/b  =  0.75.  Note  that  in  these  figures 
the  vorticity  and  streamlines  were  computed  from  phase- averaged  velocity  fields.  From 
the  figure  it  is  clear  that  the  maximum  Cm  occurred  during  predominantly  attached  flow. 
This  is  exemplified  by  Figure  20a  where  </>  =  185°  which  shows  that  the  flow  is  effectively 
attached  across  the  entire  wing  chord.  Stepping  forward  in  phase  to  Figure  20b,  the  phase- 
averaged  fields  begin  to  show  trailing  edge  initiated  separation  which  is  consistent  with  the 
stalling  behavior  of  a  relatively  thick  airfoil  profile  such  as  the  NACA  0018  utilized  here. 
By  4>  =  216°,  shown  in  Figure  20c,  the  flowfield  shows  significant  re-circulation  near  the 
trailing  edge,  which  reaches  the  leading  edge  by  <^>  =  243°.  While  the  flowfield  streamlines 
began  to  resemble  a  dynamic  stall  (leading  edge)  vortex  when  </>  =  243  this  resemblance  was 
not  indicative  of  such  a  vortex  existing  for  two  reasons.  First,  the  separation  was  trailing 
edge  initiated.  In  classical  discussions  of  dynamic  stall  the  vortex  forms  from  leading  edge 
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initiated  separation  Larsen  et  al.  [2007].  This  is  often  argued  to  be  the  result  of  a  leading 
edge  separation  bubble  bursting  event,  however  in  this  case,  no  leading  edge  separation 
bubble  forms.  Second  potential  flow  arguments  can  be  made  to  explain  why  a  coherent 
leading  edge  vortex  does  not  form  by  noting  a  particular  aspect  of  the  wing  deformation. 
Specifically,  as  noted  in  Culler  et  al.  [2017],  the  wing  deforms  in  a  spanwise  varying  fashion 
largely  defined  by  the  first  eigenmode  for  a  cantilevered  beam.  This  means  that  at  the  wing 
root,  the  airfoil  is  at  the  geometric  or  base  angle  of  attack  of  a  =  10°,  where  stall  and  large 
scale  separation  does  not  occur  in  static  tests.  Therefore,  for  a  leading  edge  vortex  to  exist 
near  the  75%  span  location,  the  vortex  filament  would  be  required  to  terminate  within  the 
fluid,  a  clear  violation  of  Helmholtz  second  theorem  or  on  the  wing  body.  Clearly  to  satisfy 
Helmholtz  theorem,  any  formed  vortices  have  to  bend  into  a  U-shape  thereby  terminating 
at  the  wing  surface. 

The  argument  that  no  dynamic  stall  vortex  is  present  is  further  supported  by  the  in¬ 
stantaneous  fields.  Specifically,  in  Figure  21,  streamlines  are  again  plotted,  however  they 
are  overlaid  on  contours  of  the  Q-criterion  which  is  discussed  by  [Hunt  et  al.,  1988]  and  is 
determined  from  Q  =  1/2(|H2|  —  \S2\)  [Jeong  and  Hussain,  1995].  In  the  plots  red  contours 
indicate  regions  of  positive  Q  and  therefore  high  rotation  while  blue  contours  indicate  neg¬ 
ative  Q  and  therefore  high  shear.  Note  that  the  figure  shows  a  sequence  of  “representative” 
instantaneous  flowfields  from  the  previously  presented  phase-averaged  sets.  Consistent  with 
the  phase-averaged  results  of  Figure  20  the  streamlines  show  a  region  of  re-circulation  be¬ 
ginning  at  the  trailing.  However,  when  examining  the  resulting  Q  fields,  it  becomes  clear 
that  the  flow  consists  of  small  scale  structure  consistent  with  a  separated  shear-layer,  rather 
than  a  coherent  dynamic  stall  vortex.  Note,  in  these  fields  Qrms  is  approximately  1.5(106). 

When  the  remaining  instantaneous  fields  are  analyzed  and  the  size  and  location  of  the 
concentrations  of  Q  are  compiled  (not  shown  here),  it  becomes  even  more  apparent  that  the 
only  small  structures  are  observed  randomly  distributed  along  a  shear  layer  emanating  from 
the  suction  surface  of  the  airfoil.  In  fact  there  is  no  coherently  identifiable  dynamic  stall 
vortex  in  any  of  the  phases.  As  a  result,  the  wing  stall  can  be  associated  with  trailing-edge 
initiated  shear  layer  separation  rather  than  the  formation  of  a  large  structure  as  is  typically 
observed  in  rigid  body  models.  Furthermore,  the  lack  of  visible  periodicity  in  the  spatial 
locations  of  the  small  scale  structures  implies  that  the  shedding  of  the  these  structures 
within  the  shear  layer  is  not  coupled  or  locked  to  the  the  frequency  of  the  LCO.  As  will  be 
discussed  in  Section  4.1.9  this  trailing-edge  initiated  shear  layer  separation  is  concentrated 
in  only  the  outer  2/3  span  of  the  flexible  wing. 

4.1.9  Spatial  Variation  of  the  Flow  Field 

As  was  shown  in  Section  4.1.8,  the  flow  separation  is  defined  by  trailing-edge  initiated  shear 
layer  separation  rather  than  the  formation  of  a  coherent,  large  scale  structure  initiating 
leading  separation.  In  order  to  understand  how  this  varies  along  the  span,  the  flowfield 
was  also  studied  at  span  locations  from  0.58  <  z/b  <  0.875  and  at  0  =  204°,  214°,  and 
224°.  Figure  22  presents  the  phase-averaged  flow  field  for  along  the  span  for  a  phase 
angle  of  (j)  =  224°,  where  streamlines  are  plotted  overlaid  on  contours  of  vorticity.  The 
figure  shows  that  the  extent  of  flow  separation  along  the  wing  chord  is  dependent  upon 
spanwise  location.  At  the  locations  z/b  =  0.58  and  0.875  as  shown  by  Figures  22a  and22f, 
respectively,  the  extent  of  flow  separation  is  significantly  reduced  compared  to  z/b  =  0.75. 
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Figure  20:  Streamlines  overlaid  on  contours  of  normalize  vorticity,  for  the  phase- 

averaged  flowhelds  around  flexible  wing  at  z/b  =  0.75  for  phase  angles  from  cj)  =  185°  to 
=  243°. 
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Figure  21:  Streamlines  overlaid  on  contours  of  Q  for  instantaneous  flowfields  around  the 
flexible  wing  at  z/b  =  0.75  for  phase  angles  from  <j)  =  185°  to  <j)  =  243°. 
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Furthermore,  at  z/b  =  0.875,  the  flow  appears  to  be  fully  attached.  This  tendency  towards 
flow  reattachnrent  is  argued  to  be  the  result  of  two  effects  depending  upon  the  spanwise 
location. 

At  inboard  spanwise  positions  the  linear  variation  in  the  wing  twist  angle  reduces  the 
effective  angle  of  attack  thus  reducing  the  severity  of  flow  separation.  Specifically,  the  wing 
root  is  fixed  to  an  angle  of  attack  of  a  =  10°  which  corresponds  to  statically  attached 
flow.  Because  of  the  twist  distribution  noted  in  Culler  et  al.  [2017],  span  locations  closer  to 
the  root  will  also  be  at  statically  attached  angles  of  attack.  Using  the  first  eigenmode  to 
estimate  this  twist  angle  at  the  spanwise  location  z/b  =  0.58  for  the  phase  angle  4>  =  224° 
it  is  found  that  the  angle  of  attack  is  expected  to  be  a  =  20.15°  based  upon  the  tip  twist 
noted  in  Figure  19.  Using  the  surface  reflection  measuring  technique  from  the  SPIV  data, 
the  angle  of  attack  at  this  location  was  measured  to  be  a  =  20.4°.  This  angle  of  attack  is 
near  the  static  stall  angle  for  a  NACA  0018  wing  [Goett  and  Bullivant,  1938]  which  would 
be  expected  to  show  some  level  of  separation.  However  because  the  wing  is  dynamically 
pitching,  the  stall  point  is  expected  to  shift  to  higher  angles  of  attack,  later  in  the  oscillation 
cycle  resulting  from  the  boundary  layer  enhancement  effects  noted  by  Ericsson  and  Reding 
[1988]  and  therefore  the  observation  of  flow  attachment  for  this  case  is  consistent. 

At  the  outboard  spanwise  postions  the  wing  tip  vortex  adds  momentum  to  the  flow  on 
the  suction  surface  and  locally  reattaches  the  flow  near  the  tip.  Specifically,  the  wing  tip 
vortex  forms  from  the  re-balancing  of  the  pressure  gradient  on  the  bottom  (pressure)  and 
top  (suction)  sides  of  the  wing.  This  induces  the  motion  of  the  air  around  the  tip,  from 
the  pressure  side  to  the  suction  side,  while  also  entraining  flow  further  away  from  the  wing 
surface  downward  onto  the  suction  side  of  the  flexible  wing  a  localized  spanwise  extent. 
This  downard  fluid  motion  helps  to  energize  and  stabilise  the  boundary  layer  near  the  wing 
tip  maintaining  attached  flow  at  extremely  high  angles  of  attack  (i.e.  a  =  30°). 

In  order  to  more  clearly  observe  how  this  region  of  separation  grows  along  the  wing  span, 
the  location  of  separation  is  identified  in  the  phase  averaged  SPIV  images  by  finding  the 
chordwise  distance  to  the  separation  point.  This  procedure  is  applied  to  each  of  the  three 
phase  locations  (</>  =  204°,  214°,  and  224°)  and  is  presented  in  Figure  23.  In  the  figure, 
the  experimentally  determined  separation  points  are  shown  by  asterisk  data  markers  and  a 
parabola  was  fit  through  each  set  (dashed  lines)  to  highlight  the  separation  distribution. 

As  the  figure  shows,  the  region  of  separation  appears  approximately  parabolic  in  distri¬ 
bution.  Furthermore,  this  region  encompasses  a  localized  area  of  the  wing  near  the  outer  2/3 
span  location.  As  the  wing  pitches  through  the  phases,  the  parabolic  region  of  separation 
grows  in  both  spanwise  and  chordwise  extent,  however  it  remains  localized.  Specifically, 
while  the  lower  bound  of  the  separation  region  remains  relatively  fixed  at  z/b  =  0.58,  the 
upper  boundary  shifts  to  higher  spanwise  positions  as  the  wing  pitches  to  larger  angles. 
This  is  because  the  wing  tip  reaches  significantly  higher  angles  of  attack  than  the  in-board 
regions  of  the  wing.  As  the  wing  pitches  up,  the  flow  separation  associated  with  these  higher 
angles  of  attack  become  too  significant  for  the  influence  of  the  tip  vortex  to  overcome. 

4.2  Rigid  Wing  Model 

Utilizing  the  cyber-physical  system  approach  discussed  in  Section  3.3,  a  comparison  is  made 
between  the  single  frequency  sinusoidally  driven,  which  will  be  subsequently  referred  to 
simply  as  ’driven’,  and  cyber-physical  stall  flutter,  subsequently  referred  to  as  ’stall  flutter’, 
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(e)  z/b  =  0.79  (f)  z/b  =  0.875 


Figure  22:  Streamlines  overlaid  on  contours  of  normalized  vorticity  from  the  phase- averaged 
flowfields  around  the  flexible  wing  from  z/b  =  0.58  to  z/b  =  0.875  at  a  phase  angle  of 
(/  =  224°. 
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- cf,  =  204 

- (j)  =  214 

- =  224 


Figure  23:  Spatial  evolution  of  the  region  of  flow  separation  along  the  span  of  the  flexible 
wing  model  at  the  phase  angles  of  (f>  =  204°,  214°,  and  224°  and  span  locations  of  z/b  =  0.58, 
0.625,  0.71,  0.75,  0.79,  and  0.875.  The  dashed  lines  indicate  the  approximate  location  of 
U  =  0  from  a  parabolic  fit  through  the  experimental  data  points  indicated  by  colored  *. 


operating  modes  for  the  NACA  0018  rigid  wing.  In  order  to  approximately  match  system 
kinematics,  cyber-physical,  stall  flutter  cases  were  captured  first.  These  cases  were  then 
post-processed  and  analyzed  to  determine  the  mean  angle  of  attack,  a,  the  pitch  oscillation 
amplitude,  A,  and  dominate  frequency,  cj,  which  were  used  to  define  the  matching  driven 
trajectory. 

4.2.1  Pitching  Moment  Comparison  -  High  Speed  (USAFA) 

From  these  comparisons,  it  was  observed  that  the  force  production  varied  quite  significantly 
between  the  driven  and  stall  flutter  modes  of  operation.  These  deviations  were  evident  when 
comparing  the  moment  coefficients  about  the  rotation  axis  (50%  chord)  where,  as  discussed 
in  Section  3.8,  the  moment  coefficient  is  presented  in  two  ways,  namely  as  the  aerodynamic 
pitching  moment  coefficient,  Cm  and  the  measured  pitching  moment  coefficient  CmtMeas.-  A 
comparison  between  a  stall  flutter  and  matched  driven  case  is  shown  in  Figure  24  for  both 
Cm  and  Cm  Meas.  versus  angle  of  attack,  a.  In  this  figure,  the  moment  coefficient  loops  were 
traversed  clockwise.  For  the  stall  flutter  case,  the  cyber-physical  system  virtual  parameters 
were  tuned  such  that  ke  =  1.55  and  rjy  =  0.015  In  both  the  stall  flutter  and 

driven  modes  of  operation,  the  tunnel  speed  was  set  to  Uqo  =  26  m/s. 

It  was  observed  that  for  this  comparison  case,  there  were  significant  differences  between 
the  moment  response  of  the  stall  flutter  and  driven  motions.  This  was  true  for  both  the 
aerodynamic  pitching  moment  coefficient,  Cm  and  the  measured  pitching  moment  coefficient 
Cm,Meas.  ■  In  the  aerodynamic  pitching  moment  coefficient,  Cm,  versus  a  plot,  shown  in 
Figure  24,  it  was  observed  that  higher  frequencies  exist.  These  frequencies  appeared  to 
cause  extra  oscillations  around  the  moment  loop.  For  the  Cm  versus  a  loop  of  the  driven 
motion,  sharp  changes  in  curvature  occurred  near  the  maximum  and  minimum  angles  of 
attack,  however  these  sharp  changes  were  not  observed  in  the  stall  flutter  case.  That 
said  it  should  be  noted  that  the  peak  angular  amplitudes  between  the  two  cases  deviated 
significantly. 
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Figure  24:  Comparison  of  Cm  (a)  and  Cm_Meas.  (b)  for  the  stall  flutter  and  driven  cases. 
Cyber-physical  system  of  Fagley  et  al.  [2016]  utilized  with  kg  =  1.55  r/y  =  0.015 
at  Uoo  =  26  m/s. 


(b)  If  Driven 


Figure  25:  Sources  of  deviation  between  the  prescribed  driven  motion  and  stall  flutter 
kinematics:  shifted  a  (a)  and  non-sinusoidal  driven  motion  as  shown  by  phase  plot  (b). 


This  deviation  was  the  result  of  two  effects:  1)  the  matching  procedure  and  2)  the  control 
authority  limitations  of  the  system.  These  two  contributions  to  the  deviation  between  stall 
flutter  and  driven  kinematics  are  visually  explained  in  Figure  25.  For  the  first  effect,  the 
driven  motion  was  defined  to  match  mean  amplitude,  time-averaged  mean  angle  of  attack, 
and  mean  frequency.  In  this  instance,  the  stall  flutter  motion  was  not  symmetric  about  its 
mean  causing  a  bias  between  the  time-averaged  mean  and  the  peak  amplitudes  as  is  shown 
in  Figure  25a.  The  second  effect  was  a  result  of  the  controller  interplay  with  the  fluid 
dynamics.  Specifically,  while  the  controller  was  prescribing  a  sinusoidal  motion,  non-linear 
aerodynamic  forces  were  encountered.  The  time  response  of  the  controller  and  total  motor 
torque  then  dictated  how  well  the  system  could  adapt  to  these  forces  without  becoming 
non-sinusoidal.  For  the  amplitudes  and  tunnel  speeds  in  this  case,  this  interplay  caused  the 
driven  motion  to  deviate  significantly  from  a  pure  sine  wave  as  shown  by  the  phase  plot  in 
Figure  25b.  Note,  a  pure  sine  wave  would  produce  a  circular  phase  plot. 

While  there  were  several  discrepancies  in  matching  the  driven  and  stall  flutter  cases 
presented  in  Figure  24,  some  physical  insight  can  still  be  gained.  Specifically,  because 
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the  driven  motion  did  not  allow  the  wing  to  move  with  respect  to  the  externally  applied 
aerodynamic  torque,  there  was  a  decoupling  of  the  wing  position  and  the  aerodynamics. 
This  can  be  examined  by  looking  at  the  force  response  near  a  =  30°  in  Figure  24a.  In  this 
region  the  stall  flutter  and  driven  trajectory  moment  coefficients  overlap.  Traversing  the 
two  moment  coefficient  loops  along  the  pitch- up,  or  clockwise,  direction,  the  stall  flutter 
trajectory  continued  to  higher  angles  of  attack,  however  the  driven  motion  was  abruptly 
forced  to  decreasing  angles  of  attack.  The  continuation  of  the  stall  flutter  motion  to  higher 
angles  of  attack  indicated  that  the  aerodynamics  were  likely  continuing  to  apply  a  large 
force.  In  contrast,  the  driven  kinematics  moved  to  a  lower  angle  of  attack  thereby  producing 
a  significant  differential  between  the  applied  motor  torque  and  the  aerodynamic  moment 
leading  to  abrupt  changes  in  moment  loop  curvature. 

When  analyzing  CmyMeas.  presented  in  Figure  24b,  the  character  of  the  pitching  moment 
changed,  significantly.  Specifically,  while  the  aerodynamic  pitching  moment  coefficient,  Cm 
showed  higher  frequency  oscillation,  the  measured  moment  coefficient  Cm)Meas.  collapsed 
to  a  relatively  simple  elliptical  response.  This  collapse  demonstrated  that  the  structure 
introduced  a  stabilizing  effect  to  the  system.  This  effect  was  not  observed  in  the  driven 
case,  however  it  should  be  noted  that  since  the  driven  kinematics  differed  from  stall  flutter, 
direct  comparisons  are  potentially  misleading. 

4.2.2  Pitching  Moment  Comparison  -  Low  Speed  (UCB) 

Due  to  the  aforementioned  issues  of  control  at  the  higher  tunnel  speeds,  the  wind  speed  was 
reduced  and  tests  were  repeated  using  the  cyber-physical  system  and  wind  tunnel  facilities 
at  the  University  of  Colorado  Boulder.  In  the  following  discussion,  a  low  speed  case  is 
analyzed  which  exhibited  highly  sinusoidal  stall  flutter  thereby  ensuring  as  equivalent  of 
a  comparison  as  possible  to  the  driven  motion.  The  kinematics  of  the  stall  flutter  and 
driven  motion  for  this  low  speed  case  is  summarized  in  Figure  26,  where  the  resulting 
moment  coefficients,  both  Crn  and  Cm>Meas.  are  shown  in  Figures  26d  and  26e.  Note  that 
the  plots  in  Figures  26b,  26c,  26d,  and  26e  show  phase  averaged  trajectories.  For  this 
stall  flutter  case,  the  cyber-physical  system  was  tuned  such  that  k$  =  1.54  N  m/(rad), 
rjv  =  0.013  N  m/(s  rad),  and  U0 0  =  15.2  m/s. 

In  Figure  26a  the  frequency  spectra  of  the  angular  position  shows  a  dominant  peak  at 

3.2  Hz,  with  the  2nd  —  4th  harmonics  also  visible.  These  frequency  peaks  are  observed  in 
both  the  stall  flutter  (black)  and  driven  (blue)  motions,  however  the  stall  flutter  case  displays 
broader  higher  harmonic  peaks,  as  compared  with  the  driven  motion.  This  potentially 
indicates  that  the  stall  flutter  motion  retains  more  energy  at  these  higher  frequencies,  though 
the  3.2  Hz  peak  is  still  dominant.  From  these  frequency  spectra  it  is  expected  that  both 
motions  are  nearly  sinusoidal.  This  is  verified  by  overlaying  a  single  sine  wave  on  the  mean 
trajectories  of  both  the  driven  and  stall  flutter  motions  in  Figure  26b.  This  overlay  shows 
that  the  driven  motion  almost  exactly  follows  the  sine  wave,  with  some  small  deviations. 
The  stall  flutter  motion  is  also  close  to  a  sine  wave,  though  larger  deviations  persist  near 
the  peak  amplitude  locations. 

While  the  trajectory  of  the  angle  of  attack  appears  largely  sinusoidal,  higher  (harmonic) 
frequencies  are  present  as  indicated  by  the  frequency  spectra.  These  higher  frequencies  are 
important  when  analyzing  the  acceleration  profiles,  which  are  shown  in  Figure  26c  for  both 
the  driven  and  stall  flutter  motions.  Specifically,  the  higher  harmonics  cause  the  accelera- 
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tion  profiles  to  deviate  significantly  away  from  a  simple  sinusoidal  profile.  From  a  purely 
mathematical  perspective  this  makes  sense  as  the  small,  periodic  deviations  from  the  si¬ 
nusoidal  position  trajectory  (see  Figure  26b)  will  be  contained  in  the  higher  harmonics  of 
the  spectra.  The  influence  of  these  higher  frequencies  is  then  amplified  through  the  double 
differentiation  of  the  position  trajectory  to  arrive  at  the  acceleration  profile.  Because  force 
is  proportional  to  acceleration,  these  frequencies  demonstrate  the  effects  of  the  nonlinear 
aerodynamic  behavior  for  the  dynamically  pitching  wing,  and  how  it  can  significantly  in¬ 
fluence  driven  or  prescribed  motions.  This  further  implies  that  even  small  deviations  from 
the  prescribed  sinusoidal  trajectory  can  have  a  strong  influence  on  the  the  aerodynamic 
force  production.  This  also  demonstrates  the  difficulty  in  performing  precise  pitching  wing 
or  dynamic  stall  experiments,  because  if  one  is  not  careful  they  can  create  significantly 
different  force  trajectories  from  similarly  prescribed  motions. 

From  the  moment  loop  comparisons  in  Figure  26d  it  is  evident  that  while  the  stall  flutter 
and  driven  signals  show  similar  behavior,  the  driven  motion  exhibits  slightly  larger  pitching 
moments  with  differing  post-stall  behavior.  To  quantify  this  difference,  a  root  mean  square 
difference  can  be  calculated  between  the  two,  phase  averaged  pitching  moment  signals  using 
Equation  (10)  where  Cm,F  refers  to  the  flutter  pitching  moment  and  Crnjj  refers  to  the 
driven.  This  RMS  difference  effectively  provides  a  standard  deviation  of  the  differences 
between  each  pitching  moment  versus  phase  signal.  This  equation  can  then  be  similarly 
written  for  Cm,Meas.-  Performing  this  calculation,  the  root  mean  square  difference  between 
the  stall  flutter  and  driven  cases  was  found  to  be  0.19  for  Cm,  and  0.20  for  Cm,Meas.-  In 
terms  of  pitching  moment  coefficient,  a  deviation  of  approximately  0.2  is  quite  significant. 

1  N 

RMS  =  -  J2  (<w(0)  - 

\  %  =  0 

4.2.3  Planar  PIV  -  Low  Speed  (UCB) 

To  further  understand  the  deviations  in  the  aerodynamic  force  production  between  the 
driven  and  stall  flutter  motions,  planar  PIV  measurements  are  presented,  focusing  at  cycle 
locations  near  the  positive  maxima  in  Cm  and  Cm,Meas.-  In  total  six  phase  angles  are 
examined  at  (f>  =  220°,  232°,  244°,  273°,  297°,  and  309°.  Figure  27  shows  where  these  phases 
are  aligned  in  the  cycle  in  terms  of  both  kinematic  and  aerodynamic  moment  profiles.  This 
figure  presents  the  phase-averaged:  aerodynamic  pitching  moment  coefficient,  Cm,  measured 
pitching  moment  coefficient,  Cm,Mea.s. ,  and  a  versus  phase  angle  for  both  the  stall  flutter 
and  driven  motions.  The  dashed  red  vertical  lines  correspond  to  each  of  the  phase  locations 
where  the  PIV  measurements  discussed  here  were  collected.  As  the  figure  shows,  the  first 
four  phases  captured  pitch-up  kinematics  while  the  final  two  phases  captured  the  post-stall, 
pitch  down  kinematics.  Specifically,  cf)  =  220°  corresponds  to  approximately  A 4>  ~  +6°  after 
the  Crn  peak  and  A  cj)  ~  —6°  before  the  CmiMeas.  peak  in  the  driven  cycle  where  this  spacing 
is  the  result  of  the  phase  resolution  of  the  collected  PIV  data,  (j)  =  232°  corresponds  to  an 
intersection  of  the  Cm  curves  for  driven  and  stall  flutter  motions  and  the  approximate  peak 
in  Cm^Meas.  for  the  stall  flutter  motion,  (j)  =  244°  corresponds  to  an  increasing  deviation  in 
Cm  between  the  stall  flutter  and  driven  motions  and  A(p  ~  +12°  after  the  peak  stall  flutter 
Cm.  4>  =  273°  captures  a  local  minima  in  Cm  and  Cm)Meas.-  Finally,  <)>  =  297°  and  309° 
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Figure  26:  Comparisons  of  frequency  spectrum  (a),  phase  averaged  angle  of  attack  (b), 
phase  averaged  angular  acceleration  (c)  versus  (j>,  and  Cm  (d)  and  Cm>Meas.  (e)  versus  angle 
of  attack  signals  for  matched  low  speed  stall  flutter  and  driven  cases. 
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(a)  Cm  (b)  Cm  ,Meas. 


(c)  Angle  of  Attack,  a 

Figure  27:  Plots  of  phase  averaged:  Cm  (a),  CmtMeas.  (b),  and  a  (c)  versus 
with  overlays  of  the  six  presented  PIV  measurement  phase  angles,  0  =  220° 

273°,  297°,  and  309°. 

correspond  to  the  formation  and  shedding  of  a  secondary  leading  edge  vortex 
resulting  flow  fields  for  each  of  these  phase  locations  are  shown  in  Figure  28  and  29  for  the 
driven  and  stall  flutter  cases  respectively.  In  these  figures  streamlines  are  plotted  with  color 
contours  of  the  normalized  spanwise  vorticity,  (nz  c)/Uoo. 

The  differences  in  pitching  moment  near  the  the  first  maximum,  as  shown  in  Figures 
26d  and  27  are  potentially  explained  from  the  PIV  measurements  at  (j>  =  220°  and  232°, 
which  show  that  LEV  formation  occurs  earlier  in  the  driven  motion  as  compared  to  stall 
flutter.  Specifically,  as  is  shown  in  Figs.  28a  and  29a,  a  vortex  begins  to  roll  up  around  the 
leading  edge  as  the  wing  pitches  up.  By  (j)  =  232°  a  coherent  LEV  has  formed  for  the  driven 
motion,  as  is  shown  in  Figure  28b.  In  contrast,  while  the  stall  flutter  motion  does  show  a 
vortex  appearing  at  (f>  =  232°  the  vortex  does  not  appear  to  be  as  strong  or  as  coherent 
through  cf>  =  244°  as  shown  in  Figure  29c. 

As  Figure  27  shows,  the  aerodynamic  pitching  moment,  Cm,  drops  significantly  after  the 
first  maximum.  This  is  likely  the  result  of  the  dynamic  stall  vortex  advecting  further  down 
the  chord  of  the  wing,  which  reduces  the  moment  arm  between  the  low-pressure  imposed 
by  the  vortex  core  on  the  wing  surface  and  the  rotation  axis.  This  finding  agrees  with 
arguments  made  by  Onoue  and  Breuer  Onoue  and  Breuer  [2016]  for  a  pitching  flat  plate 
and  is  supported  by  Figure  28c  showing  that  while  the  LEV  grows  in  size,  its  centroid  moves 
closer  to  the  rotation  axis. 

By  (j)  =  273°,  both  the  stall  flutter  and  driven  cases  pass  their  respective  maximum 
pitching  moments.  In  the  driven  case,  a  trailing  edge  vortex  appears  to  have  shed  from  the 


phase  angle 
,  232°,  244°, 


(LEV).  The 
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airfoil,  as  is  indicated  by  the  region  of  positive  vorticty,  and  a  second  LEV  begins  to  form. 
In  the  case  of  stall  flutter  (Figure  29d),  the  trailing  edge  vortex  is  still  present,  however  a 
coherent  second  LEV  has  yet  to  fully  form. 

Following  the  initial  LEV  vortex  shedding  cycles,  a  second  LEV  forms  and  sheds  dur¬ 
ing  both  the  driven  and  stall  flutter  cycles.  This  second  LEV  is  clearly  visible  at  (f>  = 
297  and  309,  given  in  Figs.  28e  and  28f  for  the  driven  motion,  and  Figs.  29e  and  29f  for  the 
stall  flutter  motion.  This  secondary  shedding  cycle  is  responsible  for  producing  an  inflec¬ 
tion  in  the  Cm  and  Cm>Meas.  versus  <f>  trajectories  that  creates  the  second  local  maxima  in 
Figs.  27a  and  27b.  Furthermore,  this  inflection  explains  why  the  small  loop  near  maximum 
a  is  observed  in  Figure  26d.  Specifically,  while  angle  of  attack  is  beginning  to  decrease,  Cm 
momentarily  changes  direction  of  curvature  and  increases.  A  similar  inflection,  but  paired 
with  a  smaller  peak  is  observed  in  the  stall  flutter  trajectory.  This  is  consistent  with  the 
smaller  amplitude  change  in  Cm  as  well  as  the  more  gradual  change  in  angle  of  attack  shown 
in  Figure  27c.  It  should  be  noted  that  this  second  shedding  cycle  is  slightly  different  than 
the  first.  Specifically,  while  the  first  LEV  cycle  showed  significant  variation  between  the 
driven  and  stall  flutter  motions,  in  the  second  cycle  the  vortex  size  and  phase  development 
are  quite  similar.  Specifically,  the  driven  motion  does  not  appear  to  lead  in  the  formation 
of  this  second  LEV,  but  instead  both  Figs.  28e  and  29e  show  similarly  sized  vortices  with 
comparable  vorticity,  encompassing  the  full  chord  of  the  airfoil.  Additionally  by  (j>  =  309° 
both  the  driven  and  stall  flutter  motions  have  shed  their  respective  LEVs  leaving  behind  a 
second  TEV.  The  stall  flutter  case  does  appear  to  show  a  slightly  smaller  second  TEV,  as 
compared  with  the  driven  case,  however  as  was  demonstrated  in  Figure  27  at  these  phase 
locations,  the  two  motions  show  similar  pitching  moment  magnitudes. 

4.2.4  Circulation  Comparison  -  Low  Speed  (UCB) 

To  more  quantitatively  analyze  the  strengths  of  the  vortices  in  the  flow,  the  circulation 
above  the  airfoil  surface  was  calculated  through  an  integration  of  vorticity  as  provided  by 
Equation  (11). 

f 2  =  T\  f  nzdA  (ii) 

cu  oo  Ja 

A  sample  integration  region  is  shown  in  Figure  30,  where  the  dashed  red  lines  show  the 
region’s  boundaries.  This  approach  was  used  to  provide  a  semi-quanitative  comparison 
between  the  driven  and  stall  flutter  cases.  Using  this  approach,  the  circulation  was  obtained 
for  the  presented  phases  and  provided  in  the  table  shown  in  Table  3.  As  the  table  shows, 
for  the  initial  LEV  formation  captured  from  220°  <  <f>  <  244°,  the  driven  motion  creates 
significantly  more  circulation.  This  result  agrees  with  the  qualitative  arguments  presented 
previously  and  the  flow  fields  shown  in  Figs.  28  and  29.  Furthermore,  this  result  makes  sense 
from  a  conservation  of  circulation  perspective.  Specifically,  through  the  first  shedding  cycle, 
the  driven  wing  reaches  higher  angles  of  attack  sooner  in  phase  than  the  driven.  From  this 
behavior,  it  would  be  expected  that  a  larger  starting  vortex  would  form,  containing  greater 
circulation  in  the  driven  motion,  than  in  the  stall  flutter  motion.  In  the  second  shedding 
cycle  the  calculation  shows  that  circulation  magnitudes  are  much  more  comparable,  though 
the  driven  motion  still  reaches  a  higher  maxima.  This  finding  is  again  consistent  with 
previous  observations. 

Due  to  the  larger  circulation  in  the  driven  cycle,  it  was  expected  that  the  driven  motion 
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Figure  28:  Streamlines  and  contours  of  vorticity,  f lz,  at  select  phases  for  the  driven  motion. 
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(e)  <p  =  297,  Stall  Flutter 
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Figure  29:  Streamlines  and  contours  of  vorticity,  12^,  at  select  phases  for  the  stall  flutter 
motion. 
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Table  3:  Calculated  circulation,  r(102), 
for  select  phases. 
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Figure  30:  Boundary  for  circulation  cal¬ 
culation. 


Figure  31:  Comparison  of  aerodynamic  differential  work  dW  normalized  by  cycle  frequency, 
uj,  for  driven  and  stall  flutter  cases. 
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operated  in  a  higher  energy  state  than  the  stall  flutter  motion.  Since  the  two  kinematic 
trajectories  were  similar,  it  was  also  expected  that  this  higher  energy  state  was  associated 
with  a  less  stable  flowfield.  In  order  to  demonstrate  this,  the  aerodynamic  differential  work: 

dW  =  —Cmda,  (12) 

was  calculated  as  discussed  by  McCroskey  McCroskey  [1981].  With  this  definition,  negative 
dW  implies  the  fluid  is  doing  work  on  the  wing  and  the  wing  will  tend  towards  instability, 
while  positive  dW  implies  the  wing  is  doing  work  on  the  fluid  and  will  tend  towards  sta¬ 
bility  McCroskey  [1981].  The  result  of  this  calculation  is  shown  in  Figure  31.  This  figure 
shows  that  the  negative  dW  in  the  driven  case  peaks  at  an  amplitude  27%  higher  than  in 
the  stall  flutter  case.  This  demonstrates  that  the  driven  motion  is  producing  a  less  stable 
flowfield,  which  is  consistent  with  the  increased  circulation  and  faster  shedding  noted  pre¬ 
viously.  The  dW  cycle  also  shows  that  during  the  pitch  up  portion  of  the  motion,  from 
<f>  ~  85°  until  (j)  ~  270°  the  aerodynamic  damping  is  negative.  This  makes  sense  as  it 
implies  the  aerodynamic  forces  are  dominating  the  structural  restoring  forces  causing  the 
wing  to  pitch  up.  During  the  remaining  phases,  the  wing  is  pitching  down  implying  that 
the  structural  restoring  forces  are  dominating  the  aerodynamic  and  the  wing  is  doing  work 
on  the  structure.  This  is  verified  by  the  positive  value  of  dW. 

5  Conclusions 

As  part  of  this  sponsored  research  effort  the  PI  and  GRA  conducted  experimental  research 
focusing  on  elucidating  the  baseline  kinematics  and  aerodynamics  associated  with  stall 
flutter  in  both  flexible  and  rigid  wing  models.  In  both  cases,  a  position-feedback  cyber¬ 
physical  (electro- mechanical)  system  was  utilized  to  simulate  the  structural  mechanics  of 
the  wing  models.  This  cyber-physical  system  allowed  for  an  in  situ  tuning  of  the  wing’s 
structural  properties  through  manipulating  the  gains  in  the  second-order  motor  control 
system.  Additionally,  the  cyber-physical  system  was  designed  to  directly  control  only  one 
degree-of-freedom,  namely  the  wing  pitch  axis  (i.e.  rotational  axis  aligned  parallel  to  the 
span).  In  the  flexible  wing,  this  allowed  for  dynamic  variations  in  the  wing  twist  distribution 
through  the  prescription  of  the  wing  tip  twist  angle  by  the  cyber-physical  system;  as  the 
wing  root  was  maintained  at  a  geometrically  fixed  angle  of  attack.  For  the  rigid  wing, 
the  entire  wing  span  was  allowed  to  rotate.  Thus  the  cyber-physical  system  prescribed  a 
uniform  pitch  angle  (or  angle  of  attack)  for  the  entire  wing  (i.e.  rigid  body  rotation). 

Flexible  Wing  Experiments  The  flexible  wing  model  consisted  of  a  finite-span,  rect¬ 
angular  wing  section  with  a  semi-span  aspect  ratio  of  six  and  composed  of  a  NACA  0018 
profile.  For  the  flexible  wing  study,  the  model’s  torsional  stiffness,  damping  and  iner¬ 
tia  were  maintained  at  kg  =  0.47  N  m  rad~l ,  rj  =  6.09  x  ICR3  N  m  s  rad _1,  and 
J  =  4.25  x  ICR4  kg  m2  rad -1,  respectively,  while  the  wind  tunnel  was  regulated  to  a 
constant  free-stream  velocity  of  U0 0  =  26  m/s.  At  these  conditions  a  torsion  dominated 
limit-cycle  oscillation  was  observed  and  was  thoroughly  classified  as  a  stall  flutter  phenom¬ 
ena  through  detailed  measurements  of  the  wing  kinematics. 

The  flexible  wing’s  kinematic  motion  was  experimentally  measured  through  the  use  of 
a  stereo-vision  motion  capture  system.  The  surface  motion  was  then  decomposed  into  the 
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principle  components  of  torsion  and  bending,  which  through  a  modal  analysis  were  projected 
onto  the  first  three-classical  modes  for  a  clamped-free  beam  in  either  torsion  or  bending 
vibration.  From  the  modal  analysis  of  the  wing  kinematics,  it  is  apparent  that  the  model 
deforms  primarily  in  torsion  and  that  this  deformation  is  principally  composed  of  the  first 
mode  of  torsional  vibration  with  small  amplitude  contributions  from  higher  order  modes. 
Furthermore,  the  motion  was  shown  to  oscillate  at  two  frequencies,  which  qualitatively 
correlated  to  the  torsional  and  bending  natural  frequencies  of  the  model  (c vg  =  5.66  Hz  and 
t Oh  =  11.45  Hz).  From  the  eigenmode  analysis  a  two  equation  empirically  fitted  model  of 
the  wing  kinematics  was  then  developed,  and  showed  very  good  agreement  with  the  raw 
stall  flutter  data.  The  stereo-vision  results  corroborated  the  wing-tip  angular  displacement 
measurements  made  with  the  optical  encoder  in  the  cyber-physical  system,  thus  validating 
this  approach  for  experimentally  replicating  and  monitoring  stall  flutter  on  the  test  article. 
While  the  empirical  model  did  a  good  job  of  capturing  the  mean  motion  observed,  its 
predictive  capability  is  limited  by  the  independently  fitted  weights  for  each  eigenmode 
included  within  the  model.  Furthermore,  additional  work  would  be  required  to  better 
model  the  time  dependent  variation  in  frequency  and  phase,  as  well  as  the  beat  observed 
in  the  bending  signal  across  the  full  testing  period. 

In  addition,  to  the  kinematic  measurements  aerodynamic  moment  and  flowfield  mea¬ 
surements  were  also  made  on  the  flexible  wing  model.  The  internal  pitching  moment  was 
computed  based  upon  the  internal  moments  in  the  flexible  wing  model  and  the  torque 
applied  by  the  cyber-physical  system.  This  calculation  demonstrated  that  the  maximum 
pitching  moment  observed  by  the  cyber-physical  system  occurred  well  before  the  peak  an¬ 
gle  of  attack  of  approximately  a  =  30°.  The  peak  pitching  moment  appeared  to  occur  just 
before  flow  separation  was  initiated  at  the  wing  trailing  edge  and  well  before  it  had  fully 
developed.  From  the  pitching  moment,  differential  work  was  calculated  showing  that  the 
majority  of  the  pitching  cycle  was  defined  by  negative  dW  which  corresponded  to  negative 
damping  and  a  tendency  to  grow  in  amplitude.  At  the  turning  locations  dW  became  slightly 
positive  indicating  a  re-stabilization  of  the  wing  system  which  bounds  the  growth  of  the  stall 
flutter  LCO  and  may  have  been  associated  with  massive  flow  separation  at  the  maximum 
angle  of  attack  and  reattachment  at  the  minimum  angle  of  attack.  Finally,  it  should  be 
noted  that  since  the  Cm  presented  did  not  incorporate  forces  at  the  second  boundary  (i.e. 
the  wing  root),  this  measure  could  only  provide  a  partial  representation  of  the  fluidic  forces. 

To  understand  the  initial  growth  of  separation  on  the  flexible  wing  surface,  stereoscopic 
particle  image  velocimetry  was  collected  at  different  phase  angles  in  the  oscillation  cycle  and 
spanwise  locations  along  the  wing.  From  the  phase  sweep  data  it  was  found  that  the  growth 
of  separation  on  the  flexible  wing  was  largely  defined  by  the  trailing  edge  initiated  separation 
of  the  shear  layer  rather  than  the  formation  of  a  leading  edge  vortex  as  is  typically  observed 
in  rigid  body  dynamic  stall.  This  separation  event  produced  small  structures  within  the 
shear  layer  that  were  not  coupled  or  locked  to  the  primary  oscillation  frequency  of  the  wing. 
From  the  spanwise  SPIV  data  it  was  found  that  the  separation  was  initiated  from  a  small 
parabolic  region  which  grew  from  the  trailing  edge  in  chordwise  and  spanwise  extent  at  as 
the  local  angle  of  attack  increased.  That  said  the  flow  separation  remained  largely  restricted 
to  the  outer  2/3  span  region  on  the  flexible  wing.  Finally,  even  with  no  data  beyond  the 
maximum  angle  of  attack  or  during  pitch  down  motion,  the  authors  strongly  believe  that 
the  data  shows  that  no  large-scale,  coherent  dynamic  stall  vortex  was  formed  by  the  flexible 
wing  during  stall  flutter. 
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Rigid  Wing  Experiments  The  second  phase  of  the  research  effort  focused  on  under¬ 
standing  the  kinematic  and  aerodynamic  differences  between  stall  flutter  and  a  prescribed  si¬ 
nusoidal  pitching  (driven)  motion  on  a  rigid  wing  section.  Tests  using  the  identical  wing  sys¬ 
tem  were  carried  out  in  wind  tunnels  at  both  USAFA  and  UCB.  It  was  found  that  even  when 
the  kinematics  were  more  closely  matched  at  lower  freestream  velocities,  U0 0  =  15.1  m/s, 
there  were  still  noticeable  differences  in  the  pitching  moments  for  the  driven  and  stall  flutter 
cases.  It  is  argued  that  these  differences  are  related  to  higher  frequencies  present  in  the 
cycle  trajectories.  Specifically,  while  these  frequencies  appeared  small  in  the  a  trajectory, 
they  are  significant  in  the  acceleration  profile  due  to  magnification  from  differentiation. 
Furthermore,  from  Newton’s  law  ( F  =  m  a),  it  is  well  understood  that  these  deviations  in 
acceleration  imply  significant  changes  in  aerodynamic  force  production  and  subsequently 
the  aerodynamic  pitching  moment  applied. 

In  order  to  understand  how  the  flowfield  was  related  to  these  differences  in  pitching 
moment,  planar  PIV  data  was  collected  near  the  point  of  maximum  pitching  moment  in 
the  cycle.  From  the  PIV  results  it  was  found  that  these  differences  could  be  attributed  to 
the  formation  and  subsequent  shedding  of  a  leading  edge  vortex  in  both  the  driven  and 
stall  flutter  cycles.  For  both  motions,  two  vortex  shedding  events  occurred,  one  near  the 
maximum  pitching  moment  coefficient  and  a  second  at  a  later  phase  angle  near  another 
local  maxima  in  the  pitching  moment  coefficient.  For  the  driven  case,  the  initial  leading 
edge  vortex  developed  faster  and  grew  larger  than  in  the  stall  flutter  case;  attaining  a 
maximum  circulation  approximately  1.6  times  the  magnitude  of  the  maximum  for  the  stall 
flutter  case.  This  was  argued  to  be  the  result  of  the  driven  wing  having  a  faster  pitch-up 
rate,  a  as  compared  to  the  stall  flutter  case.  In  order  to  conserve  circulation  in  the  field, 
this  then  required  the  formation  of  a  larger  vortex.  Interestingly,  the  secondary  formation 
and  shedding  cycle  occurred  at  similar  phases  for  both  the  driven  and  stall  flutter  motions 
even  though  the  primary  shedding  events  occurred  at  different  phases.  Throughout  the 
cycle,  both  motions  produced  similar  pitching  moments  and  amounts  of  circulation  across 
the  suction  surface  of  the  wing,  though  the  results  for  the  driven  motion  were  measurably 
larger  than  for  stall  flutter. 

Finally,  the  differential  work,  dW,  was  calculated  to  provide  a  measure  of  aerodynamic 
damping  on  the  rigid  wing  model.  This  calculation  was  argued  to  show  that  the  driven 
motion  produced  a  less  stable  flowfield  than  the  stall  flutter  motion  due  to  a  27%  increase 
in  peak  negative  dW.  It  was  further  shown  that  during  the  pitch-up  portion  of  the  cycle,  the 
aerodynamic  forces  dominated  the  structural  forces  as  indicated  by  negative  dW.  Inversely, 
and  as  expected,  the  structural  restoring  forces  dominated  the  aerodynamic  forces  during 
pitch-down  portion  of  the  cycle. 

6  Recommendations  for  Future  Work 

While  the  current  research  effort  made  significant  progress  towards  understanding  the  fun¬ 
damental  kinematics  and  flow  physics  associated  with  stall  flutter  on  both  flexible  and  rigid 
wing  sections  there  is  still  significantly  more  that  can  be  done.  The  first  major  recommen¬ 
dation  for  future  work  is  to  better  understand  the  influence  that  flow  control  has  on  both 
the  flowfield  and  the  structural  response  of  a  wing  in  a  stall  flutter  limit  cycle  oscillation. 
The  collaborators  of  the  PI  and  GRA  at  USAFA  have  already  made  significant  progress 
towards  studying  the  influence  of  flow  control  both  through  computational  simulations  and 
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experimental  studies.  Their  results  have  showed  great  promise,  but  they  have  been  limited 
to  using  only  a  single  steady-blowing  jet;  which  required  an  air  supply  and  a  solinoid  valve 
in  the  experiments.  The  PI  believes  that  future  work  utilizing  an  array  of  zero-net-mass 
flux  devices,  such  as  synthetic  jet  actuators,  has  the  potential  to  further  improve  the  ability 
of  flow  control  to  attenuate  stall  flutter  on  the  flexible  wing. 

In  addition  to  the  change  in  actuation  type,  another  challenge  associated  with  the  current 
studies  is  that  the  stall  flutter  motion  and  impact  of  flow  control  are  being  studied  only 
after  a  stall  flutter  LCO  has  stabilized.  To  actually  mitigate  the  impact  of  stall  flutter 
in  real  systems  it  needs  to  be  attenuated  before  the  oscillations  grow  into  stable  LCOs. 
As  a  result,  the  transient  growth  of  the  LCO  needs  to  be  studied  under  deterministic  and 
physically  relevant  conditions.  As  a  result  developing  a  method  of  physically  perturbing 
the  wing  system  to  deterministic  ally  grow  into  an  LCO  would  be  highly  desirable.  Once 
this  is  devised  and  demonstrated  then  the  transient  behavior  can  be  studied  and  novel  flow 
control  strategies  could  be  devised. 

Finally,  from  the  rigid  wing  work  it  was  found,  not  surprisingly,  that  the  pitch-up 
speed  and  acceleration  critically  influence  the  force  production  and  leading  edge  vortex 
development.  Furthermore,  a  wing  undergoing  stall  flutter  was  observed  to  have  a  much 
slower  pitch-up  rate  as  compared  to  a  wing  driven  in  a  sinusoidal  fashion  at  a  similar  primary 
frequency.  This  suggests  that  future  work  needs  to  be  carried  out  to  understand  how  the 
force  production  (and  flow  structures)  change  as  the  waveform  shape  is  skewed  for  a  fixed 
periodic  frequency. 

7  Accomplishments 

Multiple  accomplishments  were  achieved  throughout  the  course  of  this  research  program. 
The  scientific  and  technical  findings  of  the  program  were  discussed  in  detail  above.  The 
current  section  serves  to  summarize  the  various  accomplishments,  such  as:  honors  and 
awards  received,  presentations  given,  and  manuscripts  that  were  (or  will  be)  published  as  a 
result  of  this  work. 

7.1  Honors  and  Awards 

At  the  beginning  of  the  the  performance  period  the  GRA,  Ethan  Culler,  working  with 
the  PI  applied  for  and  was  selected  to  receive  a  National  Defense  Science  and  Engineering 
Graduate  (NDSEG)  Fellowship.  The  fellowship  is  sponsored  by  the  Air  Force  Office  of 
Scientific  Research  (AFOSR)  and  was  announced  to  the  GRA  on  April  14,  2015.  This 
fellowship  covered  the  GRAs  stipend  (salary)  and  tuition  starting  on  September  01,  2015 
and  continuing  for  three  calendar  years  to  August  31,  2018.  As  a  result  the  cooperative 
agreement  budget  was  modified  to  reallocate  the  GRA  summer  salary  funds  towards  travel 
expenses,  materials  and  supplies  for  the  project  to  enhance  the  research  funded  under 
this  cooperative  agreement.  This  budget  modification  was  subsequently  approved  by  the 
contracting  office  at  USAFA  on  August  27,  2015.  Finally,  the  later  conclusion  of  the  GRAs 
fellowship  (relative  to  the  performance  period  of  the  current  cooperative  agreement)  will 
allow  for  the  continuation  of  the  research  and  completion  of  the  technical  works  labeled  as 
“Forthcoming”  below. 
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7.2  Publications  and  Presentations 

A  brief  listing  of  all  of  the  presentations  and  publications  that  resulted  from  this  work  is 
provided  below: 

Presentations 

1.  J.  Farnsworth.  Towards  Controlling  Stall  Flutter  on  Flexible  Finite  Span 
Wings.  Invited  Seminar.  University  of  Maryland ,  College  Park,  MD,  July  20, 

2015. 

2.  E.  Culler  and  J.  Farnsworth.  Spanwise  Variation  of  Stall  Flutter  on  a  Flex¬ 
ible  NACA  0018  Finite  Span  Wing.  Oral  Presentation.  Rocky  Mountain  Fluid 
Mechanics  Research  Symposium  2015,  Boulder,  CO,  August  4,  2015. 

3.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  J.  Seidel.  Spanwise  Variation  of 
Stall  Flutter  on  a  Flexible  NACA  0018  Finite  Span  Wing.  Oral  Presentation. 
AIAA  Rocky  Mountain  Section  2015  Annual  Technical  Symposium,  Golden,  CO, 
November  6,  2015. 

4.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  T.  McLaughlin.  Spanwise  Variation 
of  Stall  Flutter  on  a  Flexible  NACA  0018  Finite  Span  Wing.  Oral  Presenation. 
54th  AIAA  Aerospace  Sciences  Meeting,  AIAA  SciTech  Forum,  San  Diego,  CA, 
January  4-8,  2016. 

5.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  J.  Seidel.  Comparison  of  driven 
and  simulated  “free”  stall  flutter  in  a  wind  tunnel.  Oral  Presentation.  69th 
Annual  Meeting  of  the  American  Physical  Society’s  Division  of  Fluid  Dynamics, 
Portland,  OR,  November  20-22,  2016. 

6.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  T.  McLaughlin.  Driven  versus  “Free 
Flutter”  Motion  of  a  NACA  0018  Finite  Span  Rigid  Wing.  Oral  Presenation. 
35th  AIAA  Applied  Aerodynamics  Conference,  AIAA  AVIATION  Forum,  Den¬ 
ver,  CO,  June  5-9,  2017. 

7.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  J.  Seidel.  The  Influence  of  Second 
Harmonic  Phase  and  Amplitude  Variation  in  Cyclically  Pitching  Wings.  Oral 
Presentation.  70th  Annual  Meeting  of  the  American  Physical  Society’s  Division 
of  Fluid  Dynamics,  Denver,  CO,  November  19-21,  2017  (Forthcoming). 

Conference  Proceedings 

1.  C.  Fagley,  J.  Seidel,  and  T.  McLaughlin,  and  J.  Farnsworth.  Aero-Servo-Elastic 
Control  of  a  Cyber-Physical  Flexible  Wing.  In  AIAA  Paper  2016-0320,  January 

2016. 

2.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  T.  McLaughlin.  Spanwise  Variation 
of  Stall  Flutter  on  a  Flexible  NACA  0018  Finite  Span  Wing.  In  AIAA  Paper 

2016- 1554 ,  January  2016. 

3.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  T.  McLaughlin.  Driven  versus  “Free 
Flutter”  Motion  of  a  NACA  0018  Finite  Span  Rigid  Wing.  In  AIAA  Paper 

2017- 4359,  June  2017. 


45 


PI:  J.  Farnsworth 


USAFA:  Flow  Control  of  Flexible  Structures 


Refereed  Publications 

1.  E.  Culler,  C.  Fagley,  J.  Seidel,  T.  McLaughlin,  and  J.  Farnsworth.  Developing 
a  reduced  order  model  from  structural  kinematic  measurements  of  a  flexible  finite 
span  wing  in  stall  flutter.  J.  Fluids  and  Structures ,  71,  pp.  56-69,  April  2017. 

2.  E.  Culler,  C.  Fagley,  J.  Seidel,  T.  McLaughlin,  and  J.  Farnsworth.  SPIV 
Flowfield  Analysis  Around  Flexible  Finite  Span  Wing  Excited  in  Stall  Flutter 
in  the  Wind  Tunnel.  In  preparation  for  submission  to  J.  Fluids  and  Structures , 
September  2017.  (Forthcoming) 

3.  E.  Culler,  C.  Fagley,  J.  Seidel,  T.  McLaughlin,  and  J.  Farnsworth.  Driven 
versus  “Free  Flutter”  Motion  of  a  NACA  0018  Finite  Span  Rigid  Wing.  In  prepa¬ 
ration  for  submission  to  J.  Fluids  and  Structures ,  November  2017.  (Forthcoming) 
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Appendices 


The  following  appendices  represent  the  conference  and  journal  publications  resulting  from 
this  work  that  are  currently  available  to  the  public  on  the  publishers  websites.  They  have 
not  been  attached  to  this  report  due  to  size  limitations,  but  are  available  upon  request  from 
the  PI  as  needed: 

1.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  T.  McLaughlin.  Spanwise  Variation  of  Stall 
Flutter  on  a  Flexible  NACA  0018  Finite  Span  Wing.  In  AIAA  Paper  2016-1554, 
January  2016. 

2.  E.  Culler,  C.  Fagley,  J.  Seidel,  T.  McLaughlin,  and  J.  Farnsworth.  Developing  a 
reduced  order  model  from  structural  kinematic  measurements  of  a  flexible  finite  span 
wing  in  stall  flutter.  J.  Fluids  and  Structures,  71,  pp.  56-69,  April  2017. 

3.  E.  Culler,  J.  Farnsworth,  C.  Fagley,  and  T.  McLaughlin.  Driven  versus  “Free  Flutter” 
Motion  of  a  NACA  0018  Finite  Span  Rigid  Wing.  In  AIAA  Paper  2017-4359,  June 
2017. 
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Symbols,  Abbreviations,  and  Acronyms 


List  of  Symbols 


A 

AR 

b 

c 

Cm 

dW 

f 

J 

ke 

M 

Q 

Q 

Rxx 

S 

t 

T 

Tx.y.z 

Coo 

Wx,i 

X 

y 

z 


amplitude 

aspect  ratio 

wing  span,  m 

wing  chord,  m 

pitching  moment  coefficient 

differential  work,  dW  =  — Cmda 

temporal  frequency,  Hz 

rotational  inertial,  kgm2 /rad 

rotational  stiffness,  Nm/rad 

rnach  number 

dynamic  pressure,  pa 

balance  between  vorticity  and  rate  of  strain  tensors,  Q  =  1/2(|02|  —  \S2\) 

normalized  product  of  autocorrelation  of  variable  X 

wing  area,  m2 

time,  s 

period,  s 

translation  (displacement)  in  the  x,y,  or  z  directions,  m 

free-stream  velocity,  m/s 

ith  mode  weight  of  variable  X 

streamwise  position,  m 

chord  normal  (vertical)  position,  m 

spanwise  position,  m 


Greek 

a 

7 

r2 

S(X) 

A(X) 

T] 

9 

IT 

P 

T 

<t> 

UJ 

Uo 

Uh 

nz 


angle  of  attack,  ° 

angular  phase  lag/lead,  ° 

circulation  aligned  in  the  z  direction,  m2 / s 

error  or  uncertainty  of  variable  X 

increment  in  variable  X 

rotational  damping,  N ms /rad 

tip-twist  angle  or  rigid  body  rotation  angle,  ° 

mathematical  constant,  3.14159 

air  density,  kg/m 3 

applied  torque,  Nrn 

periodic  phase  angle,  ° 

angular  frequency,  rad 

torsional  natural  frequency,  Hz 

bending  natural  frequency,  Hz 

vorticity  component  in  the  z  direction,  1/s 
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Subscripts 


Aero. 

=  aerodynamically  applied 

Meas. 

=  measured 

M  ot.or 

=  motor  applied 

P 

=  physical  property 

rms 

=  root  mean  square 

Struct. 

=  structurally  applied 

Tip 

=  wing  tip  condition 

V 

=  virtual  property 

List  of  Abbreviations 

Hz  = 

cycles  per  second 

kg  = 

kilogram 

m  = 

meter 

N 

newton 

pa  = 

pascal 

rad  = 

radian 

s  = 

second 

List  of  Acronyms 

DLT 

=  Direct  Linear  Transformation 

LCO 

=  Limit  Cycle  Oscillation 

NACA 

=  National  Advisory  Committee  on  Aeronautics 

PIV 

=  Particle  Image  Velocimetry 

SLA 

=  Stereolithography 

SPIV 

=  Stereoscopic  Particle  Image  Velocimetry 

UCB 

=  University  of  Colorado  Boulder 

USAFA 

=  United  States  Air  Force  Academy 
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